diff --git a/DataFormats/Detectors/TOF/include/DataFormatsTOF/Cluster.h b/DataFormats/Detectors/TOF/include/DataFormatsTOF/Cluster.h index c45567b6c9837..d5e84e7fcfead 100644 --- a/DataFormats/Detectors/TOF/include/DataFormatsTOF/Cluster.h +++ b/DataFormats/Detectors/TOF/include/DataFormatsTOF/Cluster.h @@ -16,6 +16,8 @@ #include "ReconstructionDataFormats/BaseCluster.h" #include // for base_object +#include +#include namespace o2 { @@ -27,6 +29,9 @@ namespace tof class Cluster : public o2::BaseCluster { + static constexpr float RadiusOutOfRange = 9999; // used to check if the radius was already calculated or not + static constexpr float PhiOutOfRange = 9999; // used to check if phi was already calculated or not + public: enum { kMain = 0x3FFFF, kUpLeft = 0x40000, // 2^18, 19th bit @@ -40,14 +45,16 @@ class Cluster : public o2::BaseCluster Cluster() = default; - Cluster(std::int16_t sensid, float x, float y, float z, float sy2, float sz2, float syz, float timeRaw, float time, float tot, int L0L1latency, int deltaBC); + Cluster(std::int16_t sensid, float x, float y, float z, float sy2, float sz2, float syz, double timeRaw, double time, float tot, int L0L1latency, int deltaBC); ~Cluster() = default; - float getTimeRaw() const { return mTimeRaw; } // Cluster ToF getter - void setTimeRaw(float timeRaw) { mTimeRaw = timeRaw; } // Cluster ToF setter - float getTime() const { return mTime; } // Cluster ToF getter - void setTime(float time) { mTime = time; } // Cluster ToF setter + void setBaseData(std::int16_t sensid, float x, float y, float z, float sy2, float sz2, float syz); + + double getTimeRaw() const { return mTimeRaw; } // Cluster ToF getter + void setTimeRaw(double timeRaw) { mTimeRaw = timeRaw; } // Cluster ToF setter + double getTime() const { return mTime; } // Cluster ToF getter + void setTime(double time) { mTime = time; } // Cluster ToF setter float getTot() const { return mTot; } // Cluster Charge getter void setTot(int tot) { mTot = tot; } // Cluster ToT setter int getL0L1Latency() const { return mL0L1Latency; }; // L0L1 latency @@ -55,8 +62,27 @@ class Cluster : public o2::BaseCluster int getDeltaBC() const { return mDeltaBC; }; // deltaBC void setDeltaBC(int value) { mDeltaBC = value; }; // deltaBC //float getZ() const {return mZ;} // Cluster Z - already in the definition of the cluster - float getR() const { return mR; } // Cluster Radius - float getPhi() const { return mPhi; } // Cluster Phi + float getR() + { + if (mR == RadiusOutOfRange) { + mR = TMath::Sqrt(getX() * getX() + getY() * getY() + getZ() * getZ()); + } + return mR; + } // Cluster Radius + float getPhi() + { + if (mPhi == PhiOutOfRange) { + mPhi = TMath::ATan2(getY(), getX()); + } + return mPhi; + } // Cluster Phi + int getSector() + { + if (mSector == -1) { + mSector = (TMath::ATan2(-getY(), -getX()) + TMath::Pi()) * TMath::RadToDeg() * 0.05; + } + return mSector; + } // Cluster Sector int getContributingChannels() const { return mContributingChannels; } void setContributingChannels(int contributingChannels) { mContributingChannels = contributingChannels; } @@ -104,14 +130,15 @@ class Cluster : public o2::BaseCluster private: friend class boost::serialization::access; - float mTimeRaw; // raw TOF time // CZ: in AliRoot it is a double - float mTime; // calibrated TOF time // CZ: in AliRoot it is a double + double mTimeRaw; // raw TOF time // CZ: in AliRoot it is a double + double mTime; // calibrated TOF time // CZ: in AliRoot it is a double float mTot; // Time-Over-threshold // CZ: in AliRoot it is a double int mL0L1Latency; // L0L1 latency // CZ: is it different per cluster? Checking one ESD file, it seems that it is always the same (see: /alice/data/2017/LHC17n/000280235/pass1/17000280235019.100/AliESDs.root) int mDeltaBC; // DeltaBC --> can it be a char or short? // CZ: is it different per cluster? Checking one ESD file, it seems that it can vary (see: /alice/data/2017/LHC17n/000280235/pass1/17000280235019.100/AliESDs.root) //float mZ; //! z-coordinate // CZ: to be verified if it is the same in the BaseCluster class - float mR; //! radius - float mPhi; //! phi coordinate + float mR = RadiusOutOfRange; //! radius + float mPhi = PhiOutOfRange; //! phi coordinate + int mSector = -1; //! sector number int mContributingChannels; // index of the channels that contributed to the cluster; to be read like this: // channel & 0x3FFFF -> first 18 bits to store the main channel // channel & bit19 (0x40000) -> alsoUPLEFT @@ -126,7 +153,7 @@ class Cluster : public o2::BaseCluster ClassDefNV(Cluster, 1); }; -std::ostream& operator<<(std::ostream& os, const Cluster& c); +std::ostream& operator<<(std::ostream& os, Cluster& c); } // namespace TOF } // namespace o2 #endif diff --git a/DataFormats/Detectors/TOF/src/Cluster.cxx b/DataFormats/Detectors/TOF/src/Cluster.cxx index d0d9301463aaa..f25298a7c91ce 100644 --- a/DataFormats/Detectors/TOF/src/Cluster.cxx +++ b/DataFormats/Detectors/TOF/src/Cluster.cxx @@ -14,7 +14,6 @@ #include "DataFormatsTOF/Cluster.h" #include "FairLogger.h" -#include #include #include @@ -23,14 +22,26 @@ using namespace o2::tof; ClassImp(o2::tof::Cluster); -Cluster::Cluster(std::int16_t sensid, float x, float y, float z, float sy2, float sz2, float syz, float timeRaw, float time, float tot, int L0L1Latency, int deltaBC) : o2::BaseCluster(sensid, x, y, z, sy2, sz2, syz), mTimeRaw(timeRaw), mTime(time), mTot(tot), mL0L1Latency(L0L1Latency), mDeltaBC(deltaBC), mContributingChannels(0) +Cluster::Cluster(std::int16_t sensid, float x, float y, float z, float sy2, float sz2, float syz, double timeRaw, double time, float tot, int L0L1Latency, int deltaBC) : o2::BaseCluster(sensid, x, y, z, sy2, sz2, syz), mTimeRaw(timeRaw), mTime(time), mTot(tot), mL0L1Latency(L0L1Latency), mDeltaBC(deltaBC), mContributingChannels(0) { // caching R and phi mR = TMath::Sqrt(x * x + y * y); mPhi = TMath::ATan2(y, x); + mSector = (TMath::ATan2(-getY(), -getX()) + TMath::Pi()) * TMath::RadToDeg() * 0.05; } +//______________________________________________________________________ +void Cluster::setBaseData(std::int16_t sensid, float x, float y, float z, float sy2, float sz2, float syz) +{ + setSensorID(sensid); + setXYZ(x, y, z); + setErrors(sy2, sz2, syz); + // caching R and phi + mR = TMath::Sqrt(x * x + y * y); + mPhi = TMath::ATan2(y, x); + mSector = (TMath::ATan2(-getY(), -getX()) + TMath::Pi()) * TMath::RadToDeg() * 0.05; +} //______________________________________________________________________ int Cluster::getNumOfContributingChannels() const { @@ -63,7 +74,7 @@ int Cluster::getNumOfContributingChannels() const } //______________________________________________________________________ -std::ostream& operator<<(std::ostream& os, const Cluster& c) +std::ostream& operator<<(std::ostream& os, Cluster& c) { os << (o2::BaseCluster&)c; os << " TOF cluster: raw time = " << std::scientific << c.getTimeRaw() << ", time = " << std::scientific << c.getTime() << ", Tot = " << std::scientific << c.getTot() << ", L0L1Latency = " << c.getL0L1Latency() << ", deltaBC = " << c.getDeltaBC() << ", R = " << c.getR() << ", mPhi = " << c.getPhi() << ", ContributingChannels = " << c.getNumOfContributingChannels() << "\n"; diff --git a/DataFormats/Reconstruction/CMakeLists.txt b/DataFormats/Reconstruction/CMakeLists.txt index 004a1ccceda35..5487d923baecf 100644 --- a/DataFormats/Reconstruction/CMakeLists.txt +++ b/DataFormats/Reconstruction/CMakeLists.txt @@ -7,6 +7,7 @@ set(SRCS src/BaseCluster.cxx src/TrackTPCITS.cxx src/Vertex.cxx + src/MatchInfoTOF.cxx ) Set(HEADERS @@ -14,6 +15,7 @@ Set(HEADERS include/${MODULE_NAME}/BaseCluster.h include/${MODULE_NAME}/TrackTPCITS.h include/${MODULE_NAME}/Vertex.h + include/${MODULE_NAME}/MatchInfoTOF.h ) Set(LINKDEF src/ReconstructionDataFormatsLinkDef.h) diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/MatchInfoTOF.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/MatchInfoTOF.h new file mode 100644 index 0000000000000..0d0cfdbd10869 --- /dev/null +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/MatchInfoTOF.h @@ -0,0 +1,40 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 MatchInfoTOF.h +/// \brief Class to store the output of the matching to TOF + +#ifndef ALICEO2_MATCHINFOTOF_H +#define ALICEO2_MATCHINFOTOF_H + +namespace o2 +{ +namespace dataformats +{ +class MatchInfoTOF +{ + public: + MatchInfoTOF(int indexTOFCl, float chi2) : mTOFClIndex(indexTOFCl), mChi2(chi2){}; + MatchInfoTOF() = default; + void setTOFClIndex(int index) { mTOFClIndex = index; } + int getTOFClIndex() const { return mTOFClIndex; } + + void setChi2(int chi2) { mChi2 = chi2; } + float getChi2() const { return mChi2; } + + private: + int mTOFClIndex; // index of the TOF cluster used for the matching + float mChi2; // chi2 of the pair track-TOFcluster + + // ClassDefNV(MatchInfoTOF, 1); +}; +} +} +#endif diff --git a/DataFormats/Reconstruction/src/MatchInfoTOF.cxx b/DataFormats/Reconstruction/src/MatchInfoTOF.cxx new file mode 100644 index 0000000000000..0dabe8ae8a42e --- /dev/null +++ b/DataFormats/Reconstruction/src/MatchInfoTOF.cxx @@ -0,0 +1,18 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 MatchInfoTOF.cxx +/// \brief Class to store the output of the matching to TOF + +#include "ReconstructionDataFormats/MatchInfoTOF.h" + +using namespace o2::dataformats; + +//ClassImp(o2::dataformats::MatchInfoTOF); diff --git a/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h b/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h index a62a2da8b7d08..96248e38010b5 100644 --- a/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h +++ b/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h @@ -18,6 +18,7 @@ #pragma link C++ class o2::track::TrackParCov + ; #pragma link C++ class o2::BaseCluster < float > +; #pragma link C++ class o2::dataformats::TrackTPCITS + ; +#pragma link C++ class o2::dataformats::MatchInfoTOF + ; #pragma link C++ class o2::dataformats::Vertex < int > +; #pragma link C++ class o2::dataformats::Vertex < o2::dataformats::TimeStamp < int >> +; diff --git a/DataFormats/simulation/src/SimulationDataLinkDef.h b/DataFormats/simulation/src/SimulationDataLinkDef.h index 3313d9dddc3fd..8820d35a5f030 100644 --- a/DataFormats/simulation/src/SimulationDataLinkDef.h +++ b/DataFormats/simulation/src/SimulationDataLinkDef.h @@ -37,6 +37,7 @@ #pragma link C++ struct o2::dataformats::MCTruthHeaderElement + ; #pragma link C++ class o2::dataformats::MCTruthContainer < long > +; #pragma link C++ class o2::dataformats::MCTruthContainer < o2::MCCompLabel > +; +#pragma link C++ class std::vector < o2::dataformats::MCTruthContainer < o2::MCCompLabel >> +; #pragma link C++ class std::vector < o2::MCCompLabel > +; #pragma link C++ class std::vector < o2::dataformats::MCTruthHeaderElement > +; diff --git a/Detectors/GlobalTracking/CMakeLists.txt b/Detectors/GlobalTracking/CMakeLists.txt index a6418cbeaa703..b31d81ca164d4 100644 --- a/Detectors/GlobalTracking/CMakeLists.txt +++ b/Detectors/GlobalTracking/CMakeLists.txt @@ -4,10 +4,12 @@ O2_SETUP(NAME ${MODULE_NAME}) set(SRCS src/MatchTPCITS.cxx + src/MatchTOF.cxx ) set(HEADERS include/${MODULE_NAME}/MatchTPCITS.h + include/${MODULE_NAME}/MatchTOF.h ) set(LINKDEF src/GlobalTrackingLinkDef.h) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h new file mode 100644 index 0000000000000..80137b7651fce --- /dev/null +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h @@ -0,0 +1,260 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 MatchTOF.h +/// \brief Class to perform TOF matching to global tracks +/// \author Francesco.Noferini@cern.ch, Chiara.Zampolli@cern.ch + +#ifndef ALICEO2_GLOBTRACKING_MATCHTOF_ +#define ALICEO2_GLOBTRACKING_MATCHTOF_ + +#include +#include +#include +#include +#include +#include "ReconstructionDataFormats/Track.h" +#include "ReconstructionDataFormats/TrackTPCITS.h" +#include "ReconstructionDataFormats/MatchInfoTOF.h" +#include "CommonDataFormat/EvIndex.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include "TOFBase/Geo.h" +#include "DataFormatsTOF/Cluster.h" +#include "GlobalTracking/MatchTPCITS.h" +#include "DataFormatsTPC/TrackTPC.h" + +class TTree; + +namespace o2 +{ + +namespace dataformats +{ +template +class MCTruthContainer; +} + +namespace globaltracking +{ + +///< original track in the currently loaded TPC-ITS reco output +struct TrackLocTPCITS : public o2::track::TrackParCov { + o2::dataformats::EvIndex source; ///< track origin id + timeBracket 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); +}; + +class MatchTOF +{ + using Geo = o2::tof::Geo; + using Cluster = o2::tof::Cluster; + + public: + ///< perform matching for provided input + void run(); + + ///< perform all initializations + void init(); + + ///< set tree/chain containing tracks + void setInputTreeTracks(TTree* tree) { mInputTreeTracks = tree; } + + ///< set tree/chain containing TPC tracks + void setInputTreeTPCTracks(TTree* tree) { mTreeTPCTracks = tree; } + + ///< set tree/chain containing TOF clusters + void setInputTreeTOFClusters(TTree* tree) { mTreeTOFClusters = tree; } + + ///< set output tree to write matched tracks + void setOutputTree(TTree* tr) { mOutputTree = tr; } + + ///< set input branch names for the input from the tree + void setTrackBranchName(const std::string& nm) { mTracksBranchName = nm; } + void setTPCTrackBranchName(const std::string& nm) { mTPCTracksBranchName = nm; } + void setTPCMCTruthBranchName(const std::string& nm) { mTPCMCTruthBranchName = nm; } + void setITSMCTruthBranchName(const std::string& nm) { mITSMCTruthBranchName = nm; } + void setTOFMCTruthBranchName(const std::string& nm) { mTOFMCTruthBranchName = nm; } + void setTOFClusterBranchName(const std::string& nm) { mTOFClusterBranchName = nm; } + void setOutTOFMCTruthBranchName(const std::string& nm) { mOutTOFMCTruthBranchName = nm; } + void setOutTPCMCTruthBranchName(const std::string& nm) { mOutTPCMCTruthBranchName = nm; } + void setOutITSMCTruthBranchName(const std::string& nm) { mOutITSMCTruthBranchName = nm; } + void setOutTracksBranchName(const std::string& nm) { mOutTracksBranchName = nm; } + + ///< get input branch names for the input from the tree + const std::string& getTracksBranchName() const { return mTracksBranchName; } + const std::string& getTPCTracksBranchName() const { return mTPCTracksBranchName; } + const std::string& getTPCMCTruthBranchName() const { return mTPCMCTruthBranchName; } + const std::string& getITSMCTruthBranchName() const { return mITSMCTruthBranchName; } + const std::string& getTOFMCTruthBranchName() const { return mTOFMCTruthBranchName; } + const std::string& getTOFClusterBranchName() const { return mTOFClusterBranchName; } + const std::string& getOutTOFMCTruthBranchName() const { return mOutTOFMCTruthBranchName; } + const std::string& getOutTPCMCTruthBranchName() const { return mOutTPCMCTruthBranchName; } + const std::string& getOutITSMCTruthBranchName() const { return mOutITSMCTruthBranchName; } + + ///< print settings + void print() const; + void printCandidatesTOF() const; + + ///< set time tolerance on track-TOF times comparison + void setTimeTolerance(float val) { mTimeTolerance = val; } + ///< get tolerance on track-TOF times comparison + float getTimeTolerance() const { return mTimeTolerance; } + + ///< set space tolerance on track-TOF times comparison // this in the old AliRoot was the TOF matching window + void setSpaceTolerance(float val) { mSpaceTolerance = val; } + ///< get tolerance on track-TOF 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; } + +#ifdef _ALLOW_DEBUG_TREES_ + 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); + + ///< set the name of output debug file + void setDebugTreeFileName(std::string name) + { + if (!name.empty()) { + mDebugTreeFileName = name; + } + } + + ///< get the name of output debug file + const std::string& getDebugTreeFileName() const { return mDebugTreeFileName; } + + ///< fill matching debug tree + void fillTOFmatchTree(const char* tname, int cacheTOF, int sectTOF, int plateTOF, int stripTOF, int padXTOF, int padZTOF, int cacheeTrk, int crossedStrip, int sectPropagation, int platePropagation, int stripPropagation, int padXPropagation, int padZPropagation, float resX, float resZ, float res, o2::dataformats::TrackTPCITS& trk); + void fillTOFmatchTreeWithLabels(const char* tname, int cacheTOF, int sectTOF, int plateTOF, int stripTOF, int padXTOF, int padZTOF, int cacheeTrk, int crossedStrip, int sectPropagation, int platePropagation, int stripPropagation, int padXPropagation, int padZPropagation, float resX, float resZ, float res, o2::dataformats::TrackTPCITS& trk, int TPClabelTrackID, int TPClabelEventID, int TPClabelSourceID, int ITSlabelTrackID, int ITSlabelEventID, int ITSlabelSourceID, int TOFlabelTrackID0, int TOFlabelEventID0, int TOFlabelSourceID0, int TOFlabelTrackID1, int TOFlabelEventID1, int TOFlabelSourceID1, int TOFlabelTrackID2, int TOFlabelEventID2, int TOFlabelSourceID2); + void dumpWinnerMatches(); +#endif + + private: + void attachInputTrees(); + bool prepareTracks(); + bool prepareTOFClusters(); + bool loadTracksNextChunk(); + bool loadTOFClustersNextChunk(); + + void doMatching(int sec); + void selectBestMatches(); + bool propagateToRefX(o2::track::TrackParCov& trc, float xRef /*in cm*/, float stepInCm /*in cm*/); + bool propagateToRefXWithoutCov(o2::track::TrackParCov& trc, float xRef /*in cm*/, float stepInCm /*in cm*/, float bz); + + //================================================================ + + // Data members + + bool mInitDone = false; ///< flag init already done + + float mXRef = Geo::RMIN; ///< reference radius to propage tracks for matching + + int mCurrTracksTreeEntry = -1; ///< current tracks tree entry loaded to memory + int mCurrTOFClustersTreeEntry = -1; ///< current TOF clusters tree entry loaded to memory + + bool mMCTruthON = false; ///< flag availability of MC truth + + ///========== Parameters to be set externally, e.g. from CCDB ==================== + + // to be done later + + float mTimeTolerance = 1e3; ///>>------ these are input arrays which should not be modified by the matching code + // since this info is provided by external device + std::vector* mTracksArrayInp = nullptr; ///< input tracks + std::vector* mTPCTracksArrayInp = nullptr; ///< input TPC tracks + std::vector* mTOFClustersArrayInp = nullptr; ///< input TOF clusters + + o2::dataformats::MCTruthContainer* mTOFClusLabels = nullptr; ///< input TOF clusters MC labels + std::vector mTracksLblWork; ///* mTPCLabels = nullptr; ///< TPC label of input tracks + std::vector* mITSLabels = nullptr; ///< ITS label of input tracks + + /// <<<----- + + /// mTracksWork; /// mTOFClusWork; ///, o2::constants::math::NSectors> mTracksSectIndexCache; + ///< per sector indices of TOF cluster entry in mTOFClusWork + std::array, o2::constants::math::NSectors> mTOFClusSectIndexCache; + + ///> mMatchedTracksPairs; + + /// mMatchedTracks; + std::vector> mMatchedTracks; // this is the output of the matching + std::vector mOutTOFLabels; ///< TOF label of matched tracks + std::vector mOutTPCLabels; ///< TPC label of matched tracks + std::vector mOutITSLabels; ///< ITS label of matched tracks + + int mNumOfTracks; // number of tracks to be matched + std::vector mMatchedTracksIndex; // vector of indexes of the tracks to be matched + int mNumOfClusters; // number of clusters to be matched + int* mMatchedClustersIndex = nullptr; //[mNumOfClusters] + + std::string mTracksBranchName = "TPCITS"; ///< name of branch containing input matched tracks + std::string mTPCTracksBranchName = "Tracks"; ///< name of branch containing actual TPC tracks + std::string mTPCMCTruthBranchName = "MatchTPCMCTruth"; ///< name of branch containing TPC labels + std::string mITSMCTruthBranchName = "MatchITSMCTruth"; ///< name of branch containing ITS labels + std::string mTOFMCTruthBranchName = "TOFClusterMCTruth"; ///< name of branch containing TOF clusters labels + std::string mTOFClusterBranchName = "TOFCluster"; ///< name of branch containing input ITS clusters + std::string mOutTracksBranchName = "TOFMatchInfo"; ///< name of branch containing output matched tracks + std::string mOutTOFMCTruthBranchName = "MatchTOFMCTruth"; ///< name of branch containing TOF labels for output matched tracks + std::string mOutTPCMCTruthBranchName = "MatchTPCMCTruth"; ///< name of branch containing TOF labels for output matched tracks + std::string mOutITSMCTruthBranchName = "MatchITSMCTruth"; ///< name of branch containing TOF labels for output matched tracks + +#ifdef _ALLOW_DEBUG_TREES_ + std::unique_ptr mDBGOut; + UInt_t mDBGFlags = 0; + std::string mDebugTreeFileName = "dbg_matchTOF.root"; ///< name for the debug tree file +#endif + + ///----------- aux stuff --------------/// + static constexpr float MAXSNP = 0.85; // max snp of ITS or TPC track at xRef to be matched + + TStopwatch mTimerTot; + TStopwatch mTimerDBG; + ClassDefNV(MatchTOF, 1); +}; +} // namespace globaltracking +} // namespace o2 + +#endif diff --git a/Detectors/GlobalTracking/src/GlobalTrackingLinkDef.h b/Detectors/GlobalTracking/src/GlobalTrackingLinkDef.h index 12a81052221e1..e1404038712dd 100644 --- a/Detectors/GlobalTracking/src/GlobalTrackingLinkDef.h +++ b/Detectors/GlobalTracking/src/GlobalTrackingLinkDef.h @@ -15,8 +15,15 @@ #pragma link off all functions; #pragma link C++ class o2::globaltracking::MatchTPCITS + ; +#pragma link C++ class o2::globaltracking::MatchTOF + ; #pragma link C++ class o2::globaltracking::timeBracket + ; #pragma link C++ class o2::globaltracking::TrackLocTPC + ; #pragma link C++ class o2::globaltracking::TrackLocITS + ; +#pragma link C++ class std::pair < int, o2::dataformats::MatchInfoTOF > +; +#pragma link C++ class std::vector < std::pair < int, o2::dataformats::MatchInfoTOF >> +; +#pragma link C++ class std::vector < o2::dataformats::TrackTPCITS > +; +#pragma link C++ class std::vector < o2::TPC::TrackTPC > +; +#pragma link C++ class std::vector < o2::ITS::TrackITS > +; +#pragma link C++ class std::vector < o2::tof::Cluster > +; #endif diff --git a/Detectors/GlobalTracking/src/MatchTOF.cxx b/Detectors/GlobalTracking/src/MatchTOF.cxx new file mode 100644 index 0000000000000..757b0b5555650 --- /dev/null +++ b/Detectors/GlobalTracking/src/MatchTOF.cxx @@ -0,0 +1,828 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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/Cartesian3D.h" +#include "MathUtils/Utils.h" +#include "CommonConstants/MathConstants.h" +#include "CommonConstants/PhysicsConstants.h" +#include "DetectorsBase/GeometryManager.h" + +#include +#include +#include +#include +#include "DataFormatsParameters/GRPObject.h" + +#include "GlobalTracking/MatchTOF.h" + +using namespace o2::globaltracking; +using timeEst = o2::dataformats::TimeStampWithError; +using evIdx = o2::dataformats::EvIndex; + +ClassImp(MatchTOF); + +//______________________________________________ +void MatchTOF::run() +{ + ///< running the matching + + if (!mInitDone) { + LOG(FATAL) << "init() was not done yet"; + } + + mTimerTot.Start(); + + // we load all TOF clusters (to be checked if we need to split per time frame) + prepareTOFClusters(); + + // we do the matching per entry of the TPCITS matched tracks tree + while (mCurrTracksTreeEntry + 1 < mInputTreeTracks->GetEntries()) { // we add "+1" because mCurrTracksTreeEntry starts from -1, and it is incremented in loadTracksNextChunk which is called by prepareTracks + LOG(DEBUG) << "Number of entries in track tree = " << mCurrTracksTreeEntry; + mMatchedTracks.clear(); + mOutTOFLabels.clear(); + mOutTPCLabels.clear(); + mOutITSLabels.clear(); + prepareTracks(); + + /* Uncomment for local debug + Printf("*************** Printing the tracks before starting the matching"); + + // printing the tracks + std::array globalPosTmp; + int totTracks = 0; + + for (int sec = o2::constants::math::NSectors; sec--;) { + Printf("\nsector %d", sec); + auto& cacheTrkTmp = mTracksSectIndexCache[sec]; // array of cached tracks indices for this sector; reminder: they are ordered in time! + for (int itrk = 0; itrk < cacheTrkTmp.size(); itrk++){ + auto& trc = mTracksWork[cacheTrkTmp[itrk]]; + trc.getXYZGlo(globalPosTmp); + LOG(INFO) << "Track" << totTracks << " [in this sector it is the " << itrk << "]: Global coordinates After propagating to 371 cm: globalPos[0] = " << globalPosTmp[0] << ", globalPos[1] = " << globalPosTmp[1] << ", globalPos[2] = " << globalPosTmp[2]; + LOG(INFO) << "The phi angle is " << TMath::ATan2(globalPosTmp[1], globalPosTmp[0]); + totTracks++; + } + } + */ + + for (int sec = o2::constants::math::NSectors; sec--;) { + LOG(INFO) << "Doing matching for sector " << sec << "..."; + doMatching(sec); + LOG(INFO) << "...done. Now check the best matches"; + selectBestMatches(); + } + if (0) { // enabling this creates very verbose output + mTimerTot.Stop(); + printCandidatesTOF(); + mTimerTot.Start(false); + } + mOutputTree->Fill(); + } + +#ifdef _ALLOW_DEBUG_TREES_ + mDBGOut.reset(); +#endif + + mTimerTot.Stop(); + printf("Timing:\n"); + printf("Total: "); + mTimerTot.Print(); +} + +//______________________________________________ +void MatchTOF::init() +{ + ///< initizalizations + + if (mInitDone) { + LOG(ERROR) << "Initialization was already done"; + return; + } + + attachInputTrees(); + + // create output branch + if (mOutputTree) { + mOutputTree->Branch(mOutTracksBranchName.data(), &mMatchedTracks); + LOG(INFO) << "Matched tracks will be stored in " << mOutTracksBranchName << " branch of tree " + << mOutputTree->GetName(); + if (mMCTruthON) { + mOutputTree->Branch(mOutTPCMCTruthBranchName.data(), &mOutITSLabels); + LOG(INFO) << "ITS Tracks Labels branch: " << mOutITSMCTruthBranchName; + mOutputTree->Branch(mOutITSMCTruthBranchName.data(), &mOutTPCLabels); + LOG(INFO) << "TPC Tracks Labels branch: " << mOutTPCMCTruthBranchName; + mOutputTree->Branch(mOutTOFMCTruthBranchName.data(), &mOutTOFLabels); + LOG(INFO) << "TOF Tracks Labels branch: " << mOutTOFMCTruthBranchName; + } + + } else { + LOG(ERROR) << "Output tree is not attached, matched tracks will not be stored"; + } + +#ifdef _ALLOW_DEBUG_TREES_ + // debug streamer + if (mDBGFlags) { + mDBGOut = std::make_unique(mDebugTreeFileName.data(), "recreate"); + } +#endif + + mInitDone = true; + + { + mTimerTot.Stop(); + mTimerTot.Reset(); + } + + print(); +} + +//______________________________________________ +void MatchTOF::print() const +{ + ///< print the settings + + LOG(INFO) << "****** component for the matching of tracks to TOF clusters ******"; + if (!mInitDone) { + LOG(INFO) << "init is not done yet - nothing to print"; + return; + } + + LOG(INFO) << "MC truth: " << (mMCTruthON ? "on" : "off"); + LOG(INFO) << "Time tolerance: " << mTimeTolerance; + LOG(INFO) << "Space tolerance: " << mSpaceTolerance; + LOG(INFO) << "SigmaTimeCut: " << mSigmaTimeCut; + + LOG(INFO) << "**********************************************************************"; +} + +//______________________________________________ +void MatchTOF::printCandidatesTOF() const +{ + ///< print the candidates for the matching +} + +//______________________________________________ +void MatchTOF::attachInputTrees() +{ + ///< attaching the input tree + + if (!mInputTreeTracks) { + LOG(FATAL) << "Input tree with tracks is not set"; + } + + if (!mTreeTPCTracks) { + LOG(FATAL) << "TPC tracks data input tree is not set"; + } + + if (!mTreeTOFClusters) { + LOG(FATAL) << "TOF clusters data input tree is not set"; + } + + // input tracks (this is the pais of ITS-TPC matches) + + if (!mInputTreeTracks->GetBranch(mTracksBranchName.data())) { + LOG(FATAL) << "Did not find tracks branch " << mTracksBranchName << " in the input tree"; + } + mInputTreeTracks->SetBranchAddress(mTracksBranchName.data(), &mTracksArrayInp); + LOG(INFO) << "Attached tracks " << mTracksBranchName << " branch with " << mInputTreeTracks->GetEntries() + << " entries"; + + // input TOF clusters + + if (!mTreeTOFClusters->GetBranch(mTOFClusterBranchName.data())) { + LOG(FATAL) << "Did not find TOF clusters branch " << mTOFClusterBranchName << " in the input tree"; + } + mTreeTOFClusters->SetBranchAddress(mTOFClusterBranchName.data(), &mTOFClustersArrayInp); + LOG(INFO) << "Attached TOF clusters " << mTOFClusterBranchName << " branch with " << mTreeTOFClusters->GetEntries() + << " entries"; + + // is there MC info available ? + if (mTreeTOFClusters->GetBranch(mTOFMCTruthBranchName.data())) { + mTreeTOFClusters->SetBranchAddress(mTOFMCTruthBranchName.data(), &mTOFClusLabels); + LOG(INFO) << "Found TOF Clusters MCLabels branch " << mTOFMCTruthBranchName; + } + if (mInputTreeTracks->GetBranch(mTPCMCTruthBranchName.data())) { + mInputTreeTracks->SetBranchAddress(mTPCMCTruthBranchName.data(), &mTPCLabels); + LOG(INFO) << "Found TPC tracks MCLabels branch " << mTPCMCTruthBranchName.data(); + } + if (mInputTreeTracks->GetBranch(mITSMCTruthBranchName.data())) { + mInputTreeTracks->SetBranchAddress(mITSMCTruthBranchName.data(), &mITSLabels); + LOG(INFO) << "Found ITS tracks MCLabels branch " << mITSMCTruthBranchName.data(); + } + + mMCTruthON = (mTOFClusLabels && mTPCLabels && mITSLabels); + mCurrTracksTreeEntry = -1; + mCurrTOFClustersTreeEntry = -1; +} + +//______________________________________________ +bool MatchTOF::prepareTracks() +{ + ///< prepare the tracks that we want to match to TOF + + if (!loadTracksNextChunk()) { + return false; + } + + mNumOfTracks = mTracksArrayInp->size(); + if (mNumOfTracks == 0) + return false; // no tracks to be matched + mMatchedTracksIndex.resize(mNumOfTracks); + std::fill(mMatchedTracksIndex.begin(), mMatchedTracksIndex.end(), -1); // initializing all to -1 + + // copy the track params, propagate to reference X and build sector tables + mTracksWork.clear(); + mTracksWork.reserve(mNumOfTracks); + if (mMCTruthON) { + mTracksLblWork.clear(); + mTracksLblWork.reserve(mNumOfTracks); + } + for (int sec = o2::constants::math::NSectors; sec--;) { + mTracksSectIndexCache[sec].clear(); + mTracksSectIndexCache[sec].reserve(100 + 1.2 * mNumOfTracks / o2::constants::math::NSectors); + } + + // getting Bz (mag field) + auto o2field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + float bzField = o2field->solenoidField(); // magnetic field in kGauss + + Printf("\n\nWe have %d tracks to try to match to TOF", mNumOfTracks); + int nNotPropagatedToTOF = 0; + for (int it = 0; it < mNumOfTracks; it++) { + o2::dataformats::TrackTPCITS& trcOrig = (*mTracksArrayInp)[it]; // TODO: check if we cannot directly use the o2::track::TrackParCov class instead of o2::dataformats::TrackTPCITS, and then avoid the casting below; this is the track at the vertex + std::array globalPos; + + // create working copy of track param + mTracksWork.emplace_back(trcOrig); //, mCurrTracksTreeEntry, it); + // make a copy of the TPC track that we have to propagate + //o2::TPC::TrackTPC* trc = new o2::TPC::TrackTPC(trcTPCOrig); // this would take the TPCout track + auto& trc = mTracksWork.back(); // with this we take the TPCITS track propagated to the vertex + + // propagate to matching Xref + trc.getXYZGlo(globalPos); + LOG(DEBUG) << "Global coordinates Before propagating to 371 cm: globalPos[0] = " << globalPos[0] << ", globalPos[1] = " << globalPos[1] << ", globalPos[2] = " << globalPos[2]; + LOG(DEBUG) << "Radius xy Before propagating to 371 cm = " << TMath::Sqrt(globalPos[0] * globalPos[0] + globalPos[1] * globalPos[1]); + LOG(DEBUG) << "Radius xyz Before propagating to 371 cm = " << TMath::Sqrt(globalPos[0] * globalPos[0] + globalPos[1] * globalPos[1] + globalPos[2] * globalPos[2]); + if (!propagateToRefXWithoutCov(trc, mXRef, 2, bzField)) { // we first propagate to 371 cm without considering the covariance matrix + nNotPropagatedToTOF++; + continue; + } + + // the "rough" propagation worked; now we can propagate considering also the cov matrix + if (!propagateToRefX(trc, mXRef, 2) || TMath::Abs(trc.getZ()) > Geo::MAXHZTOF) { // we check that the propagation with the cov matrix worked; CHECK: can it happen that it does not if the progataion without the errors succeeded? + nNotPropagatedToTOF++; + continue; + } + + trc.getXYZGlo(globalPos); + + LOG(DEBUG) << "Global coordinates After propagating to 371 cm: globalPos[0] = " << globalPos[0] << ", globalPos[1] = " << globalPos[1] << ", globalPos[2] = " << globalPos[2]; + LOG(DEBUG) << "Radius xy After propagating to 371 cm = " << TMath::Sqrt(globalPos[0] * globalPos[0] + globalPos[1] * globalPos[1]); + LOG(DEBUG) << "Radius xyz After propagating to 371 cm = " << TMath::Sqrt(globalPos[0] * globalPos[0] + globalPos[1] * globalPos[1] + globalPos[2] * globalPos[2]); + LOG(DEBUG) << "The track will go to sector " << o2::utils::Angle2Sector(TMath::ATan2(globalPos[1], globalPos[0])); + + mTracksSectIndexCache[o2::utils::Angle2Sector(TMath::ATan2(globalPos[1], globalPos[0]))].push_back(it); + //delete trc; // Check: is this needed? + } + + LOG(INFO) << "Total number of tracks = " << mNumOfTracks << ", Number of tracks that failed to be propagated to TOF = " << nNotPropagatedToTOF; + + // sort tracks in each sector according to their time (increasing in time) + for (int sec = o2::constants::math::NSectors; sec--;) { + auto& indexCache = mTracksSectIndexCache[sec]; + LOG(INFO) << "Sorting sector" << sec << " | " << indexCache.size() << " tracks"; + if (!indexCache.size()) + continue; + std::sort(indexCache.begin(), indexCache.end(), [this](int a, int b) { + auto& trcA = mTracksWork[a]; + auto& trcB = mTracksWork[b]; + return ((trcA.getTimeMUS().getTimeStamp() - mSigmaTimeCut * trcA.getTimeMUS().getTimeStampError()) - (trcB.getTimeMUS().getTimeStamp() - mSigmaTimeCut * trcB.getTimeMUS().getTimeStampError()) < 0.); + }); + } // loop over tracks of single sector + + // Uncomment for local debug + /* + // printing the tracks + std::array globalPos; + int itmp = 0; + for (int sec = o2::constants::math::NSectors; sec--;) { + Printf("sector %d", sec); + auto& cacheTrk = mTracksSectIndexCache[sec]; // array of cached tracks indices for this sector; reminder: they are ordered in time! + for (int itrk = 0; itrk < cacheTrk.size(); itrk++){ + itmp++; + auto& trc = mTracksWork[cacheTrk[itrk]]; + trc.getXYZGlo(globalPos); + printf("Track %d: Global coordinates After propagating to 371 cm: globalPos[0] = %f, globalPos[1] = %f, globalPos[2] = %f\n", itrk, globalPos[0], globalPos[1], globalPos[2]); + // Printf("The phi angle is %f", TMath::ATan2(globalPos[1], globalPos[0])); + } + } + Printf("we have %d tracks",itmp); + */ + + return true; +} +//______________________________________________ +bool MatchTOF::prepareTOFClusters() +{ + ///< prepare the tracks that we want to match to TOF + + // copy the track params, propagate to reference X and build sector tables + mTOFClusWork.clear(); + // mTOFClusWork.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--;) { + mTOFClusSectIndexCache[sec].clear(); + //mTOFClusSectIndexCache[sec].reserve(100 + 1.2 * mNumOfClusters / o2::constants::math::NSectors); // we cannot do this, we don't have mNumOfClusters yet + } + + mNumOfClusters = 0; + while (loadTOFClustersNextChunk()) { + int nClusterInCurrentChunk = mTOFClustersArrayInp->size(); + LOG(DEBUG) << "nClusterInCurrentChunk = " << nClusterInCurrentChunk; + mNumOfClusters += nClusterInCurrentChunk; + for (int it = 0; it < nClusterInCurrentChunk; it++) { + Cluster& clOrig = (*mTOFClustersArrayInp)[it]; + + // create working copy of track param + mTOFClusWork.emplace_back(clOrig); + auto& cl = mTOFClusWork.back(); + // cache work track index + mTOFClusSectIndexCache[o2::utils::Angle2Sector(cl.getPhi())].push_back(mTOFClusWork.size() - 1); + } + } + + // sort clusters in each sector according to their time (increasing in time) + for (int sec = o2::constants::math::NSectors; sec--;) { + auto& indexCache = mTOFClusSectIndexCache[sec]; + LOG(INFO) << "Sorting sector" << sec << " | " << indexCache.size() << " TOF clusters"; + if (!indexCache.size()) + continue; + std::sort(indexCache.begin(), indexCache.end(), [this](int a, int b) { + auto& clA = mTOFClusWork[a]; + auto& clB = mTOFClusWork[b]; + return (clA.getTime() - clB.getTime()) < 0.; + }); + } // loop over TOF clusters of single sector + + if (mMatchedClustersIndex) + delete[] mMatchedClustersIndex; + mMatchedClustersIndex = new int[mNumOfClusters]; + std::fill_n(mMatchedClustersIndex, mNumOfClusters, -1); // initializing all to -1 + + return true; +} + +//_____________________________________________________ +bool MatchTOF::loadTracksNextChunk() +{ + ///< load next chunk of tracks to be matched to TOF + while (++mCurrTracksTreeEntry < mInputTreeTracks->GetEntries()) { + mInputTreeTracks->GetEntry(mCurrTracksTreeEntry); + LOG(INFO) << "Loading tracks entry " << mCurrTracksTreeEntry << " -> " << mTracksArrayInp->size() + << " tracks"; + if (!mTracksArrayInp->size()) { + continue; + } + return true; + } + --mCurrTracksTreeEntry; + return false; +} +//______________________________________________ +bool MatchTOF::loadTOFClustersNextChunk() +{ + ///< load next chunk of clusters to be matched to TOF + printf("Loading TOF clusters: number of entries in tree = %lld\n", mTreeTOFClusters->GetEntries()); + while (++mCurrTOFClustersTreeEntry < mTreeTOFClusters->GetEntries()) { + mTreeTOFClusters->GetEntry(mCurrTOFClustersTreeEntry); + LOG(DEBUG) << "Loading TOF clusters entry " << mCurrTOFClustersTreeEntry << " -> " << mTOFClustersArrayInp->size() + << " clusters"; + LOG(INFO) << "Loading TOF clusters entry " << mCurrTOFClustersTreeEntry << " -> " << mTOFClustersArrayInp->size() + << " clusters"; + if (!mTOFClustersArrayInp->size()) { + continue; + } + return true; + } + --mCurrTOFClustersTreeEntry; + return false; +} +//______________________________________________ +void MatchTOF::doMatching(int sec) +{ + ///< do the real matching per sector + + mMatchedTracksPairs.clear(); // new sector + + //uncomment for local debug + /* + // printing the tracks + std::array globalPosTmp; + Printf("sector %d", sec); + auto& cacheTrkTmp = mTracksSectIndexCache[sec]; // array of cached tracks indices for this sector; reminder: they are ordered in time! + for (int itrk = 0; itrk < cacheTrkTmp.size(); itrk++){ + auto& trc = mTracksWork[cacheTrkTmp[itrk]]; + trc.getXYZGlo(globalPosTmp); + printf("Track %d: Global coordinates After propagating to 371 cm: globalPos[0] = %f, globalPos[1] = %f, globalPos[2] = %f\n", itrk, globalPosTmp[0], globalPosTmp[1], globalPosTmp[2]); + Printf("The phi angle is %f", TMath::ATan2(globalPosTmp[1], globalPosTmp[0])); + } + */ + auto& cacheTOF = mTOFClusSectIndexCache[sec]; // array of cached TOF cluster indices for this sector; reminder: they are ordered in time! + auto& cacheTrk = mTracksSectIndexCache[sec]; // array of cached tracks indices for this sector; reminder: they are ordered in time! + int nTracks = cacheTrk.size(), nTOFCls = cacheTOF.size(); + LOG(INFO) << "Matching sector " << sec << ": number of tracks: " << nTracks << ", number of TOF clusters: " << nTOFCls; + if (!nTracks || !nTOFCls) { + return; + } + int itof0 = 0; // starting index in TOF 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 + 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]; + + 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& trefTrk = mTracksWork[cacheTrk[itrk]]; + float minTrkTime = (trefTrk.getTimeMUS().getTimeStamp() - mSigmaTimeCut * trefTrk.getTimeMUS().getTimeStampError()) * 1.E6; // minimum time in ps + float maxTrkTime = (trefTrk.getTimeMUS().getTimeStamp() + mSigmaTimeCut * trefTrk.getTimeMUS().getTimeStampError()) * 1.E6; // maximum time in ps + int istep = 1; // number of steps + float step = 0.1; // 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])); + */ + +#ifdef _ALLOW_DEBUG_TREES_ + (*mDBGOut) << "propOK" + << "track=" << trefTrk << "\n"; +#endif + + // 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 }; // TOF detector id at the current propagation point + while (propagateToRefX(trefTrk, mXRef + istep * step, step) && nStripsCrossedInPropagation <= 2 && mXRef + istep * step < Geo::RMAX) { + if (0 && istep % 100 == 0) { + printf("istep = %d, currentPosition = %f \n", istep, mXRef + istep * step); + } + 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; + Geo::getPadDxDyDz(posFloat, detIdTemp, deltaPosTemp); + + // uncomment below only for local debug; this will produce A LOT of output - one print per propagation step + //Printf("detIdTemp[0] = %d, detIdTemp[1] = %d, detIdTemp[2] = %d, detIdTemp[3] = %d, detIdTemp[4] = %d", detIdTemp[0], detIdTemp[1], detIdTemp[2], detIdTemp[3], detIdTemp[4]); + if (detIdTemp[2] != -1 && nStripsCrossedInPropagation == 0) { // print in case you have a useful propagation + LOG(DEBUG) << "*********** We have crossed a strip during propagation!*********"; + LOG(DEBUG) << "Global coordinates: pos[0] = " << pos[0] << ", pos[1] = " << pos[1] << ", pos[2] = " << pos[2]; + LOG(DEBUG) << "detIdTemp[0] = " << detIdTemp[0] << ", detIdTemp[1] = " << detIdTemp[1] << ", detIdTemp[2] = " << detIdTemp[2] << ", detIdTemp[3] = " << detIdTemp[3] << ", detIdTemp[4] = " << detIdTemp[4]; + LOG(DEBUG) << "deltaPosTemp[0] = " << deltaPosTemp[0] << ", deltaPosTemp[1] = " << deltaPosTemp[1] << " deltaPosTemp[2] = " << deltaPosTemp[2]; + } else { + LOG(DEBUG) << "*********** We have NOT crossed a strip during propagation!*********"; + LOG(DEBUG) << "Global coordinates: pos[0] = " << pos[0] << ", pos[1] = " << pos[1] << ", pos[2] = " << pos[2]; + LOG(DEBUG) << "detIdTemp[0] = " << detIdTemp[0] << ", detIdTemp[1] = " << detIdTemp[1] << ", detIdTemp[2] = " << detIdTemp[2] << ", detIdTemp[3] = " << detIdTemp[3] << ", detIdTemp[4] = " << detIdTemp[4]; + LOG(DEBUG) << "deltaPosTemp[0] = " << deltaPosTemp[0] << ", deltaPosTemp[1] = " << deltaPosTemp[1] << " deltaPosTemp[2] = " << deltaPosTemp[2]; + } + istep++; + // check if after the propagation we are in a TOF strip + if (detIdTemp[2] != -1) { // we ended in a TOF strip + LOG(DEBUG) << "nStripsCrossedInPropagation = " << nStripsCrossedInPropagation << ", detId[nStripsCrossedInPropagation][0] = " << detId[nStripsCrossedInPropagation][0] << ", detIdTemp[0] = " << detIdTemp[0] << ", detId[nStripsCrossedInPropagation][1] = " << detId[nStripsCrossedInPropagation][1] << ", detIdTemp[1] = " << detIdTemp[1] << ", detId[nStripsCrossedInPropagation][2] = " << detId[nStripsCrossedInPropagation][2] << ", detIdTemp[2] = " << detIdTemp[2]; + 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++; + } + //Printf("nStepsInsideSameStrip[nStripsCrossedInPropagation-1] = %d", nStepsInsideSameStrip[nStripsCrossedInPropagation-1]); + if (nStepsInsideSameStrip[nStripsCrossedInPropagation - 1] == 0) { + detId[nStripsCrossedInPropagation - 1][0] = detIdTemp[0]; + detId[nStripsCrossedInPropagation - 1][1] = detIdTemp[1]; + detId[nStripsCrossedInPropagation - 1][2] = detIdTemp[2]; + detId[nStripsCrossedInPropagation - 1][3] = detIdTemp[3]; + detId[nStripsCrossedInPropagation - 1][4] = detIdTemp[4]; + deltaPos[nStripsCrossedInPropagation - 1][0] = deltaPosTemp[0]; + deltaPos[nStripsCrossedInPropagation - 1][1] = deltaPosTemp[1]; + deltaPos[nStripsCrossedInPropagation - 1][2] = deltaPosTemp[2]; + nStepsInsideSameStrip[nStripsCrossedInPropagation - 1]++; + } else { // a further propagation step in the same strip -> update info (we sum up on all matching with strip - we will divide for the number of steps a bit below) + deltaPos[nStripsCrossedInPropagation - 1][0] += deltaPosTemp[0] + (detIdTemp[4] - detId[nStripsCrossedInPropagation - 1][4]) * Geo::XPAD; // residual in x + deltaPos[nStripsCrossedInPropagation - 1][1] += deltaPosTemp[1]; // residual in y + deltaPos[nStripsCrossedInPropagation - 1][2] += deltaPosTemp[2] + (detIdTemp[3] - detId[nStripsCrossedInPropagation - 1][3]) * Geo::ZPAD; // residual in z + nStepsInsideSameStrip[nStripsCrossedInPropagation - 1]++; + } + } + } + LOG(DEBUG) << "while done, we propagated track " << itrk << " in %d strips" << nStripsCrossedInPropagation; + + // uncomment for debug purposes, to check tracks that did not cross any strip + /* + if (nStripsCrossedInPropagation == 0) { + auto labelTPCNoStripsCrossed = mTPCLabels->at(mTracksSectIndexCache[sec][itrk]); + Printf("The current track (index = %d) never crossed a strip", cacheTrk[itrk]); + Printf("TrackID = %d, EventID = %d, SourceID = %d", labelTPCNoStripsCrossed.getTrackID(), labelTPCNoStripsCrossed.getEventID(), labelTPCNoStripsCrossed.getSourceID()); + printf("Global coordinates: pos[0] = %f, pos[1] = %f, pos[2] = %f\n", pos[0], pos[1], pos[2]); + printf("detIdTemp[0] = %d, detIdTemp[1] = %d, detIdTemp[2] = %d, detIdTemp[3] = %d, detIdTemp[4] = %d\n", detIdTemp[0], detIdTemp[1], detIdTemp[2], detIdTemp[3], detIdTemp[4]); + printf("deltaPosTemp[0] = %f, deltaPosTemp[1] = %f, deltaPosTemp[2] = %f\n", deltaPosTemp[0], deltaPosTemp[1], deltaPosTemp[2]); + } + */ + + 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]; + LOG(DEBUG) << "matched strip " << imatch << ": deltaPos[0] = " << deltaPos[imatch][0] << ", deltaPos[1] = " << deltaPos[imatch][1] << ", deltaPos[2] = " << deltaPos[imatch][2] << ", residual (x, z) = " << TMath::Sqrt(deltaPos[imatch][0] * deltaPos[imatch][0] + deltaPos[imatch][2] * deltaPos[imatch][2]); + } + + if (nStripsCrossedInPropagation == 0) { + continue; // the track never hit a TOF strip during the propagation + } + Printf("We will check now the %d TOF clusters", nTOFCls); + bool foundCluster = false; + auto labelTPC = (*mTPCLabels)[mTracksSectIndexCache[sec][itrk]]; + for (auto itof = itof0; itof < nTOFCls; itof++) { + // printf("itof = %d\n", itof); + auto& trefTOF = mTOFClusWork[cacheTOF[itof]]; + // 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); + /* This part is commented out for now, as we don't want to have any check on the time enabled + 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"); + itof0 = itof+1; // but for the next track that we will check, we will ignore this cluster (the time is anyway too small) + //continue; + } + if (trefTOF.getTime() > maxTrkTime) { // no more TOF clusters can be matched to this track + // break; + } + */ + int mainChannel = trefTOF.getMainContributingChannel(); + int indices[5]; + Geo::getVolumeIndices(mainChannel, indices); + const auto& labelsTOF = mTOFClusLabels->getLabels(mTOFClusSectIndexCache[indices[0]][itof]); + int trackIdTOF; + int eventIdTOF; + int sourceIdTOF; + for (auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation; iPropagation++) { + LOG(DEBUG) << "TOF Cluster [" << itof << ", " << cacheTOF[itof] << "]: indices = " << indices[0] << ", " << indices[1] << ", " << indices[2] << ", " << indices[3] << ", " << indices[4]; + LOG(DEBUG) << "Propagated Track [" << itrk << ", " << cacheTrk[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; // 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; // 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; +#ifdef _ALLOW_DEBUG_TREES_ + fillTOFmatchTree("match0", cacheTOF[itof], indices[0], indices[1], indices[2], indices[3], indices[4], cacheTrk[itrk], iPropagation, detId[iPropagation][0], detId[iPropagation][1], detId[iPropagation][2], detId[iPropagation][3], detId[iPropagation][4], resX, resZ, res, trefTrk); +#endif + int tofLabelTrackID[3] = { -1, -1, -1 }; + int tofLabelEventID[3] = { -1, -1, -1 }; + int tofLabelSourceID[3] = { -1, -1, -1 }; + for (int ilabel = 0; ilabel < labelsTOF.size(); ilabel++) { + tofLabelTrackID[ilabel] = labelsTOF[ilabel].getTrackID(); + tofLabelEventID[ilabel] = labelsTOF[ilabel].getEventID(); + tofLabelSourceID[ilabel] = labelsTOF[ilabel].getSourceID(); + } + //auto labelTPC = mTPCLabels->at(mTracksSectIndexCache[indices[0]][itrk]); + auto labelITS = (*mITSLabels)[mTracksSectIndexCache[indices[0]][itrk]]; + fillTOFmatchTreeWithLabels("matchPossibleWithLabels", cacheTOF[itof], indices[0], indices[1], indices[2], indices[3], indices[4], cacheTrk[itrk], iPropagation, detId[iPropagation][0], detId[iPropagation][1], detId[iPropagation][2], detId[iPropagation][3], detId[iPropagation][4], resX, resZ, res, trefTrk, labelTPC.getTrackID(), labelTPC.getEventID(), labelTPC.getSourceID(), labelITS.getTrackID(), labelITS.getEventID(), labelITS.getSourceID(), tofLabelTrackID[0], tofLabelEventID[0], tofLabelSourceID[0], tofLabelTrackID[1], tofLabelEventID[1], tofLabelSourceID[1], tofLabelTrackID[2], tofLabelEventID[2], tofLabelSourceID[2]); + 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! + fillTOFmatchTree("match1", cacheTOF[itof], indices[0], indices[1], indices[2], indices[3], indices[4], cacheTrk[itrk], iPropagation, detId[iPropagation][0], detId[iPropagation][1], detId[iPropagation][2], detId[iPropagation][3], detId[iPropagation][4], resX, resZ, res, trefTrk); + + fillTOFmatchTreeWithLabels("matchOkWithLabels", cacheTOF[itof], indices[0], indices[1], indices[2], indices[3], indices[4], cacheTrk[itrk], iPropagation, detId[iPropagation][0], detId[iPropagation][1], detId[iPropagation][2], detId[iPropagation][3], detId[iPropagation][4], resX, resZ, res, trefTrk, labelTPC.getTrackID(), labelTPC.getEventID(), labelTPC.getSourceID(), labelITS.getTrackID(), labelITS.getEventID(), labelITS.getSourceID(), tofLabelTrackID[0], tofLabelEventID[0], tofLabelSourceID[0], tofLabelTrackID[1], tofLabelEventID[1], tofLabelSourceID[1], tofLabelTrackID[2], tofLabelEventID[2], tofLabelSourceID[2]); + + if (res < mSpaceTolerance) { // matching ok! + LOG(DEBUG) << "MATCHING FOUND: We have a match! between track " << mTracksSectIndexCache[indices[0]][itrk] << " and TOF cluster " << mTOFClusSectIndexCache[indices[0]][itof]; + foundCluster = true; + mMatchedTracksPairs.emplace_back(std::make_pair(mTracksSectIndexCache[indices[0]][itrk], o2::dataformats::MatchInfoTOF(mTOFClusSectIndexCache[indices[0]][itof], chi2))); // TODO: check if this is correct! + for (int ilabel = 0; ilabel < labelsTOF.size(); ilabel++) { + LOG(DEBUG) << "TOF label " << ilabel << ": trackID = " << labelsTOF[ilabel].getTrackID() << ", eventID = " << labelsTOF[ilabel].getEventID() << ", sourceID = " << labelsTOF[ilabel].getSourceID(); + } + LOG(DEBUG) << "TPC label of the track: trackID = " << labelTPC.getTrackID() << ", eventID = " << labelTPC.getEventID() << ", sourceID = " << labelTPC.getSourceID(); + LOG(DEBUG) << "ITS label of the track: trackID = " << labelITS.getTrackID() << ", eventID = " << labelITS.getEventID() << ", sourceID = " << labelITS.getSourceID(); + fillTOFmatchTreeWithLabels("matchOkWithLabelsInSpaceTolerance", cacheTOF[itof], indices[0], indices[1], indices[2], indices[3], indices[4], cacheTrk[itrk], iPropagation, detId[iPropagation][0], detId[iPropagation][1], detId[iPropagation][2], detId[iPropagation][3], detId[iPropagation][4], resX, resZ, res, trefTrk, labelTPC.getTrackID(), labelTPC.getEventID(), labelTPC.getSourceID(), labelITS.getTrackID(), labelITS.getEventID(), labelITS.getSourceID(), tofLabelTrackID[0], tofLabelEventID[0], tofLabelSourceID[0], tofLabelTrackID[1], tofLabelEventID[1], tofLabelSourceID[1], tofLabelTrackID[2], tofLabelEventID[2], tofLabelSourceID[2]); + } + } + } + if (!foundCluster) + LOG(DEBUG) << "We did not find any TOF cluster for track " << cacheTrk[itrk] << " (label = " << labelTPC.getTrackID() << ", pt = " << trefTrk.getPt(); + } + return; +} + +//______________________________________________ +void MatchTOF::selectBestMatches() +{ + ///< define the track-TOFcluster pair per sector + + // first, we sort according to the chi2 + std::sort(mMatchedTracksPairs.begin(), mMatchedTracksPairs.end(), [this](std::pair a, std::pair b) { return (a.second.getChi2() < b.second.getChi2()); }); + int i = 0; + // then we take discard the pairs if their track or cluster was already matched (since they are ordered in chi2, we will take the best matching) + for (const std::pair& matchingPair : mMatchedTracksPairs) { + if (mMatchedTracksIndex[matchingPair.first] != -1) { // the track was already filled + continue; + } + if (mMatchedClustersIndex[matchingPair.second.getTOFClIndex()] != -1) { // the track was already filled + continue; + } + mMatchedTracksIndex[matchingPair.first] = mMatchedTracks.size(); // index of the MatchInfoTOF correspoding to this track + mMatchedClustersIndex[matchingPair.second.getTOFClIndex()] = mMatchedTracksIndex[matchingPair.first]; // index of the track that was matched to this cluster + mMatchedTracks.push_back(matchingPair); // array of MatchInfoTOF + const auto& labelTPC = (*mTPCLabels)[matchingPair.first]; + LOG(DEBUG) << "labelTPC: trackID = " << labelTPC.getTrackID() << ", eventID = " << labelTPC.getEventID() << ", sourceID = " << labelTPC.getSourceID(); + const auto& labelITS = (*mITSLabels)[matchingPair.first]; + LOG(DEBUG) << "labelITS: trackID = " << labelITS.getTrackID() << ", eventID = " << labelITS.getEventID() << ", sourceID = " << labelITS.getSourceID(); + const auto& labelsTOF = mTOFClusLabels->getLabels(matchingPair.second.getTOFClIndex()); + bool labelOk = false; // whether we have found or not the same TPC label of the track among the labels of the TOF cluster + for (int ilabel = 0; ilabel < labelsTOF.size(); ilabel++) { + LOG(DEBUG) << "TOF label " << ilabel << ": trackID = " << labelsTOF[ilabel].getTrackID() << ", eventID = " << labelsTOF[ilabel].getEventID() << ", sourceID = " << labelsTOF[ilabel].getSourceID(); + if (labelsTOF[ilabel].getTrackID() == labelTPC.getTrackID() && labelsTOF[ilabel].getEventID() == labelTPC.getEventID() && labelsTOF[ilabel].getSourceID() == labelTPC.getSourceID() && !labelOk) { // if we find one TOF cluster label that is the same as the TPC one, we are happy - even if it is not the first one + mOutTOFLabels.push_back(labelsTOF[ilabel]); + labelOk = true; + } + } + if (!labelOk) { + // we have not found the track label among those associated to the TOF cluster --> fake match! We will associate the label of the main channel, but negative + o2::MCCompLabel fakeTOFlabel; + fakeTOFlabel.set(-labelsTOF[0].getTrackID(), labelsTOF[0].getEventID(), labelsTOF[0].getSourceID()); + mOutTOFLabels.push_back(fakeTOFlabel); + } + mOutTPCLabels.push_back(labelTPC); + mOutITSLabels.push_back(labelITS); + i++; + } +} + +//______________________________________________ +bool MatchTOF::propagateToRefX(o2::track::TrackParCov& trc, float xRef, float stepInCm) +{ + // propagate track to matching reference X + 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, o2::constants::physics::MassPionCharged, MAXSNP, stepInCm, 0.); + 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::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, o2::constants::physics::MassPionCharged, MAXSNP, stepInCm, 0.); + } + // 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 MatchTOF::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::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 +} + +#ifdef _ALLOW_DEBUG_TREES_ +//______________________________________________ +void MatchTOF::setDebugFlag(UInt_t flag, bool on) +{ + ///< set debug stream flag + if (on) { + mDBGFlags |= flag; + } else { + mDBGFlags &= ~flag; + } +} + +//_________________________________________________________ +void MatchTOF::fillTOFmatchTree(const char* trname, int cacheTOF, int sectTOF, int plateTOF, int stripTOF, int padXTOF, int padZTOF, int cacheeTrk, int crossedStrip, int sectPropagation, int platePropagation, int stripPropagation, int padXPropagation, int padZPropagation, float resX, float resZ, float res, o2::dataformats::TrackTPCITS& trk) +{ + ///< fill debug tree for TOF tracks matching check + + mTimerDBG.Start(false); + + // Printf("************** Filling the debug tree with %d, %d, %d, %d, %d, %d, %d, %d, %d, %d, %d, %d, %d, %f, %f, %f", cacheTOF, sectTOF, plateTOF, stripTOF, padXTOF, padZTOF, cacheeTrk, crossedStrip, sectPropagation, platePropagation, stripPropagation, padXPropagation, padZPropagation, resX, resZ, res); + + (*mDBGOut) << trname + << "clusterTOF=" << cacheTOF << "sectTOF=" << sectTOF << "plateTOF=" << plateTOF << "stripTOF=" << stripTOF << "padXTOF=" << padXTOF << "padZTOF=" << padZTOF + << "crossedStrip=" << crossedStrip << "sectPropagation=" << sectPropagation << "platePropagation=" << platePropagation << "stripPropagation=" << stripPropagation << "padXPropagation=" << padXPropagation + << "resX=" << resX << "resZ=" << resZ << "res=" << res << "track=" << trk << "\n"; + mTimerDBG.Stop(); +} + +//_________________________________________________________ +void MatchTOF::fillTOFmatchTreeWithLabels(const char* trname, int cacheTOF, int sectTOF, int plateTOF, int stripTOF, int padXTOF, int padZTOF, int cacheeTrk, int crossedStrip, int sectPropagation, int platePropagation, int stripPropagation, int padXPropagation, int padZPropagation, float resX, float resZ, float res, o2::dataformats::TrackTPCITS& trk, int TPClabelTrackID, int TPClabelEventID, int TPClabelSourceID, int ITSlabelTrackID, int ITSlabelEventID, int ITSlabelSourceID, int TOFlabelTrackID0, int TOFlabelEventID0, int TOFlabelSourceID0, int TOFlabelTrackID1, int TOFlabelEventID1, int TOFlabelSourceID1, int TOFlabelTrackID2, int TOFlabelEventID2, int TOFlabelSourceID2) +{ + ///< fill debug tree for TOF tracks matching check + + mTimerDBG.Start(false); + + (*mDBGOut) << trname + << "clusterTOF=" << cacheTOF << "sectTOF=" << sectTOF << "plateTOF=" << plateTOF << "stripTOF=" << stripTOF << "padXTOF=" << padXTOF << "padZTOF=" << padZTOF + << "crossedStrip=" << crossedStrip << "sectPropagation=" << sectPropagation << "platePropagation=" << platePropagation << "stripPropagation=" << stripPropagation << "padXPropagation=" << padXPropagation + << "resX=" << resX << "resZ=" << resZ << "res=" << res << "track=" << trk + << "TPClabelTrackID=" << TPClabelTrackID << "TPClabelEventID=" << TPClabelEventID << "TPClabelSourceID=" << TPClabelSourceID + << "ITSlabelTrackID=" << ITSlabelTrackID << "ITSlabelEventID=" << ITSlabelEventID << "ITSlabelSourceID=" << ITSlabelSourceID + << "TOFlabelTrackID0=" << TOFlabelTrackID0 << "TOFlabelEventID0=" << TOFlabelEventID0 << "TOFlabelSourceID0=" << TOFlabelSourceID0 + << "TOFlabelTrackID1=" << TOFlabelTrackID1 << "TOFlabelEventID1=" << TOFlabelEventID1 << "TOFlabelSourceID1=" << TOFlabelSourceID1 + << "TOFlabelTrackID2=" << TOFlabelTrackID2 << "TOFlabelEventID2=" << TOFlabelEventID2 << "TOFlabelSourceID2=" << TOFlabelSourceID2 + << "\n"; + mTimerDBG.Stop(); +} +#endif diff --git a/Detectors/TOF/base/include/TOFBase/Digit.h b/Detectors/TOF/base/include/TOFBase/Digit.h index 2e8b9edba5721..e56320c82a117 100644 --- a/Detectors/TOF/base/include/TOFBase/Digit.h +++ b/Detectors/TOF/base/include/TOFBase/Digit.h @@ -11,7 +11,6 @@ #ifndef ALICEO2_TOF_DIGIT_H_ #define ALICEO2_TOF_DIGIT_H_ -#include #include #include "Rtypes.h" @@ -21,12 +20,12 @@ namespace o2 { namespace tof { /// \class Digit /// \brief TOF digit implementation -class Digit : public o2::dataformats::TimeStamp +class Digit { public: Digit() = default; - Digit(Double_t time, Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t label = -1); + Digit(Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t label = -1); ~Digit() = default; /// Get global ordering key made of @@ -52,7 +51,7 @@ class Digit : public o2::dataformats::TimeStamp void printStream(std::ostream &stream) const; - void merge(Double_t time, Int_t tdc, Int_t tot); + void merge(Int_t tdc, Int_t tot); void getPhiAndEtaIndex(int& phi, int& eta); diff --git a/Detectors/TOF/base/include/TOFBase/Geo.h b/Detectors/TOF/base/include/TOFBase/Geo.h index b81297da4c8fa..26be900a21f2d 100644 --- a/Detectors/TOF/base/include/TOFBase/Geo.h +++ b/Detectors/TOF/base/include/TOFBase/Geo.h @@ -94,7 +94,8 @@ class Geo static constexpr Float_t DEADTIME = 25E+03; // Single channel dead time (ps) static constexpr Float_t DEADTIMETDC = DEADTIME/TDCBIN; ///< Single channel TDC dead time (ps) static constexpr Float_t MATCHINGWINDOW = TDCBIN * 8192; // Matching window (ps) 2^13=8192 - static constexpr Float_t READOUTWINDOW = 1000; // Readout window (ns) + // static constexpr Float_t READOUTWINDOW = 1000; // Readout window (ns) + static constexpr Float_t READOUTWINDOW = 1e9; // Readout window (ns) - now put 1s for DPL to work fine, but it will be 29e3 static constexpr Float_t READOUTWINDOW_INV = 1. / READOUTWINDOW; // Readout window (ns) static constexpr Float_t ANGLES[NPLATES][NMAXNSTRIP] = { // Strip Tilt Angles diff --git a/Detectors/TOF/base/src/Digit.cxx b/Detectors/TOF/base/src/Digit.cxx index 869972905cce1..4c89cb11b8c46 100644 --- a/Detectors/TOF/base/src/Digit.cxx +++ b/Detectors/TOF/base/src/Digit.cxx @@ -17,8 +17,8 @@ using namespace o2::tof; ClassImp(o2::tof::Digit); -Digit::Digit(Double_t time, Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t label) - : o2::dataformats::TimeStamp(time), mChannel(channel), mTDC(tdc), mTOT(tot), mBC(bc), mLabel(label), mIsUsedInCluster(kFALSE) +Digit::Digit(Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t label) + : mChannel(channel), mTDC(tdc), mTOT(tot), mBC(bc), mLabel(label), mIsUsedInCluster(kFALSE) { } @@ -26,8 +26,7 @@ Digit::Digit(Double_t time, Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t void Digit::printStream(std::ostream& stream) const { - stream << "TOF Digit: Channel " << mChannel << " TDC " << mTDC << " TOT " << mTOT << " Time " << getTimeStamp() - << "Bunch Crossing index" << mBC << " Label " << mLabel << "\n"; + stream << "TOF Digit: Channel " << mChannel << " TDC " << mTDC << " TOT " << mTOT << "Bunch Crossing index" << mBC << " Label " << mLabel << "\n"; } //______________________________________________________________________ @@ -40,14 +39,13 @@ std::ostream& operator<<(std::ostream& stream, const Digit& digi) //______________________________________________________________________ -void Digit::merge(Double_t time, Int_t tdc, Int_t tot) +void Digit::merge(Int_t tdc, Int_t tot) { // merging two digits if (tdc < mTDC) { mTDC = tdc; - setTimeStamp(time); // TODO: adjust TOT } else { // TODO: adjust TOT diff --git a/Detectors/TOF/base/src/Geo.cxx b/Detectors/TOF/base/src/Geo.cxx index 64763f69b19e2..af8bfe08bcff7 100644 --- a/Detectors/TOF/base/src/Geo.cxx +++ b/Detectors/TOF/base/src/Geo.cxx @@ -12,6 +12,7 @@ #include "TGeoManager.h" #include "TMath.h" #include "FairLogger.h" +#include "DetectorsBase/GeometryManager.h" ClassImp(o2::tof::Geo); @@ -104,7 +105,7 @@ void Geo::Init() void Geo::getVolumePath(const Int_t* ind, Char_t* path) { //-------------------------------------------------------------------- - // This function returns the colume path of a given pad + // This function returns the volume path of a given pad //-------------------------------------------------------------------- Int_t sector = ind[0]; @@ -157,8 +158,13 @@ void Geo::getPos(Int_t* det, Float_t* pos) Char_t path[200]; getVolumePath(det, path); if (!gGeoManager) { - LOG(ERROR) << " no TGeo!" << "\n"; + LOG(ERROR) << " no TGeo! Loading it" + << "\n"; + o2::Base::GeometryManager::loadGeometry(); } + FILE* ciccio = fopen("TOF_geo.txt", "w"); + fprintf(ciccio, "path = %s, gGeoManager = %p", path, gGeoManager); + fclose(ciccio); gGeoManager->cd(path); TGeoHMatrix global; global = *gGeoManager->GetCurrentMatrix(); @@ -321,7 +327,7 @@ Int_t Geo::getStripNumberPerSM(Int_t iplate, Int_t istrip) void Geo::fromGlobalToSector(Float_t* pos, Int_t isector) { if (isector == -1) { - LOG(ERROR) << "Sector Index not valid (-1)\n"; + //LOG(ERROR) << "Sector Index not valid (-1)\n"; return; } @@ -369,9 +375,7 @@ Int_t Geo::fromPlateToStrip(Float_t* pos, Int_t iplate) step[1] = getHeights(iplate, istrip); step[2] = -getDistances(iplate, istrip); translate(posLoc2, step); - rotateToStrip(posLoc2, iplate, istrip); - if ((TMath::Abs(posLoc2[0]) <= STRIPLENGTH * 0.5) && (TMath::Abs(posLoc2[1]) <= HSTRIPY * 0.5) && (TMath::Abs(posLoc2[2]) <= WCPCBZ * 0.5)) { step[0] = -0.5 * NPADX * XPAD; diff --git a/Detectors/TOF/base/src/TOFBaseLinkDef.h b/Detectors/TOF/base/src/TOFBaseLinkDef.h index 4ce5aa3b7b478..3cef9f5920fb7 100644 --- a/Detectors/TOF/base/src/TOFBaseLinkDef.h +++ b/Detectors/TOF/base/src/TOFBaseLinkDef.h @@ -17,5 +17,6 @@ #pragma link C++ class o2::tof::Geo+; #pragma link C++ class o2::tof::Digit+; #pragma link C++ class vector+; +#pragma link C++ class vector < vector < o2::tof::Digit >> +; #endif diff --git a/Detectors/TOF/prototyping/findLabels.C b/Detectors/TOF/prototyping/findLabels.C new file mode 100644 index 0000000000000..e1cd2439adc00 --- /dev/null +++ b/Detectors/TOF/prototyping/findLabels.C @@ -0,0 +1,85 @@ +void findLabels(int itrk, int ientry) +{ + + // macro to find the labels of a TPCITS track and the corresponding TOF cluster + + // getting the ITSTPCtracks + TFile* fmatchITSTPC = new TFile("o2match_itstpc.root"); + TTree* matchTPCITS = (TTree*)fmatchITSTPC->Get("matchTPCITS"); + std::vector* mTracksArrayInp = new std::vector; + matchTPCITS->SetBranchAddress("TPCITS", &mTracksArrayInp); + matchTPCITS->GetEntry(ientry); + + // getting the TPC tracks + TFile* ftracksTPC = new TFile("tpctracks.root"); + TTree* tpcTree = (TTree*)ftracksTPC->Get("events"); + std::vector* mTPCTracksArrayInp = new std::vector; + tpcTree->SetBranchAddress("Tracks", &mTPCTracksArrayInp); + o2::dataformats::MCTruthContainer* mcTPC = new o2::dataformats::MCTruthContainer(); + tpcTree->SetBranchAddress("TPCTracksMCTruth", &mcTPC); + + // getting the ITS tracks + TFile* ftracksITS = new TFile("o2trac_its.root"); + TTree* itsTree = (TTree*)ftracksITS->Get("o2sim"); + std::vector* mITSTracksArrayInp = new std::vector; + itsTree->SetBranchAddress("ITSTrack", &mITSTracksArrayInp); + o2::dataformats::MCTruthContainer* mcITS = new o2::dataformats::MCTruthContainer(); + itsTree->SetBranchAddress("ITSTrackMCTruth", &mcITS); + + // getting the TOF clusters + TFile* fclustersTOF = new TFile("tofclusters.root"); + TTree* tofClTree = (TTree*)fclustersTOF->Get("o2sim"); + std::vector* mTOFClustersArrayInp = new std::vector; + tofClTree->SetBranchAddress("TOFCluster", &mTOFClustersArrayInp); + o2::dataformats::MCTruthContainer* mcTOF = new o2::dataformats::MCTruthContainer(); + tofClTree->SetBranchAddress("TOFClusterMCTruth", &mcTOF); + + tpcTree->GetEntry(0); + tofClTree->GetEntry(0); + + o2::dataformats::TrackTPCITS trackITSTPC = mTracksArrayInp->at(itrk); + const o2::dataformats::EvIndex& evIdxTPC = trackITSTPC.getRefTPC(); + Printf("matched TPCtrack: eventID = %d, indexID = %d", evIdxTPC.getEvent(), evIdxTPC.getIndex()); + const o2::dataformats::EvIndex& evIdxITS = trackITSTPC.getRefITS(); + Printf("matched ITStrack: eventID = %d, indexID = %d", evIdxITS.getEvent(), evIdxITS.getIndex()); + itsTree->GetEntry(evIdxITS.getEvent()); + + // getting the TPC labels + const auto& labelsTPC = mcTPC->getLabels(evIdxTPC.getIndex()); + for (int ilabel = 0; ilabel < labelsTPC.size(); ilabel++) { + Printf("TPC label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsTPC[ilabel].getTrackID(), labelsTPC[ilabel].getEventID(), labelsTPC[ilabel].getSourceID()); + } + + // getting the ITS labels + const auto& labelsITS = mcITS->getLabels(evIdxITS.getIndex()); + for (int ilabel = 0; ilabel < labelsITS.size(); ilabel++) { + Printf("ITS label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsITS[ilabel].getTrackID(), labelsITS[ilabel].getEventID(), labelsITS[ilabel].getSourceID()); + } + + // now checking if we have a corresponding cluster + bool found = false; + for (int tofClIndex = 0; tofClIndex < mTOFClustersArrayInp->size(); tofClIndex++) { + o2::tof::Cluster tofCluster = mTOFClustersArrayInp->at(tofClIndex); + const auto& labelsTOF = mcTOF->getLabels(tofClIndex); + int trackIdTOF; + int eventIdTOF; + int sourceIdTOF; + for (int ilabel = 0; ilabel < labelsTOF.size(); ilabel++) { + if ((labelsTPC[0].getTrackID() == labelsTOF[ilabel].getTrackID() && labelsTPC[0].getEventID() == labelsTOF[ilabel].getEventID() && labelsTPC[0].getSourceID() == labelsTOF[ilabel].getSourceID()) || (labelsITS[0].getTrackID() == labelsTOF[ilabel].getTrackID() && labelsITS[0].getEventID() == labelsTOF[ilabel].getEventID() && labelsITS[0].getSourceID() == labelsTOF[ilabel].getSourceID())) { + Printf("The corresponding TOF cluster is %d", tofClIndex); + Printf("TOF label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsTOF[ilabel].getTrackID(), labelsTOF[ilabel].getEventID(), labelsTOF[ilabel].getSourceID()); + found = true; + int nContrChannels = tofCluster.getNumOfContributingChannels(); + int mainContrChannel = tofCluster.getMainContributingChannel(); + int* indices = new int[5]; + o2::tof::Geo::getVolumeIndices(mainContrChannel, indices); + Printf("main contributing channel: sector = %d, plate = %d, strip = %d, padz = %d, padx = %d", indices[0], indices[1], indices[2], indices[3], indices[4]); + break; + } + } + } + if (!found) + Printf("No TOF cluster corresponding to this track was found"); + + return; +} diff --git a/Detectors/TOF/prototyping/findTOFclusterFromLabel.C b/Detectors/TOF/prototyping/findTOFclusterFromLabel.C new file mode 100644 index 0000000000000..b5ae58b0120b8 --- /dev/null +++ b/Detectors/TOF/prototyping/findTOFclusterFromLabel.C @@ -0,0 +1,115 @@ +void findTOFclusterFromLabel(int trackID, int eventID = 0, int sourceID = 0) +{ + + // macro to find the labels of a TPCITS track and the corresponding TOF cluster + + // getting the TOF clusters + TFile* fclustersTOF = new TFile("tofclusters.root"); + TTree* tofClTree = (TTree*)fclustersTOF->Get("o2sim"); + std::vector* mTOFClustersArrayInp = new std::vector; + tofClTree->SetBranchAddress("TOFCluster", &mTOFClustersArrayInp); + o2::dataformats::MCTruthContainer* mcTOF = new o2::dataformats::MCTruthContainer(); + tofClTree->SetBranchAddress("TOFClusterMCTruth", &mcTOF); + + tofClTree->GetEntry(0); + + // now checking if we have a corresponding cluster + bool found = false; + for (int tofClIndex = 0; tofClIndex < mTOFClustersArrayInp->size(); tofClIndex++) { + o2::tof::Cluster tofCluster = mTOFClustersArrayInp->at(tofClIndex); + const auto& labelsTOF = mcTOF->getLabels(tofClIndex); + int trackIdTOF; + int eventIdTOF; + int sourceIdTOF; + for (int ilabel = 0; ilabel < labelsTOF.size(); ilabel++) { + if (trackID == labelsTOF[ilabel].getTrackID() && eventID == labelsTOF[ilabel].getEventID() && sourceID == labelsTOF[ilabel].getSourceID()) { + Printf("The corresponding TOF cluster is %d", tofClIndex); + Printf("TOF label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsTOF[ilabel].getTrackID(), labelsTOF[ilabel].getEventID(), labelsTOF[ilabel].getSourceID()); + found = true; + int nContrChannels = tofCluster.getNumOfContributingChannels(); + int mainContrChannel = tofCluster.getMainContributingChannel(); + int* indices = new int[5]; + o2::tof::Geo::getVolumeIndices(mainContrChannel, indices); + Printf("main contributing channel: sector = %d, plate = %d, strip = %d, padz = %d, padx = %d", indices[0], indices[1], indices[2], indices[3], indices[4]); + break; + } + } + } + if (!found) + Printf("No TOF cluster corresponding to this label was found"); + + TFile* fdigitsTOF = new TFile("tofdigits.root"); + TTree* tofDigTree = (TTree*)fdigitsTOF->Get("o2sim"); + std::vector>* mTOFDigitsArrayInp = nullptr; + tofDigTree->SetBranchAddress("TOFDigit", &mTOFDigitsArrayInp); + std::vector>* mcTOFDig = nullptr; + tofDigTree->SetBranchAddress("TOFDigitMCTruth", &mcTOFDig); + bool foundInDigits = false; + Printf("Now looking in the digits"); + for (int ientry = 0; ientry < tofDigTree->GetEntries(); ientry++) { + //Printf("\n\nEntry in tree %d", ientry); + tofDigTree->GetEntry(ientry); + for (int iVect = 0; iVect < mTOFDigitsArrayInp->size(); iVect++) { + //Printf("\nEntry in vector of digits and MC truth %d", iVect); + std::vector digitVector = mTOFDigitsArrayInp->at(iVect); + o2::dataformats::MCTruthContainer digitMCTruth = mcTOFDig->at(iVect); + for (int iDig = 0; iDig < digitVector.size(); iDig++) { + //Printf("Digit %d", iDig); + o2::tof::Digit digit = digitVector.at(iDig); + int digitLabel = digit.getLabel(); + gsl::span mcArray = digitMCTruth.getLabels(digitLabel); + for (int j = 0; j < static_cast(mcArray.size()); j++) { + //printf("checking element %d in the array of labels\n", j); + auto label = digitMCTruth.getElement(digitMCTruth.getMCTruthHeader(digitLabel).index + j); + //printf("TrackID = %d, EventID = %d, SourceID = %d\n", label.getTrackID(), label.getEventID(), label.getSourceID()); + if (label.getTrackID() == trackID && label.getEventID() == eventID && label.getSourceID() == sourceID) { + Printf("We found the label that we were looking for! tree entry = %d, vector entry = %d, digit = %d", ientry, iVect, iDig); + foundInDigits = true; + } + } + } + } + } + if (!foundInDigits) + Printf("The label was NEVER found in the digits"); + + TFile* fKine = new TFile("o2sim.root"); + TTree* tKine = (TTree*)fKine->Get("o2sim"); + std::vector* mcArr = nullptr; + tKine->SetBranchAddress("MCTrack", &mcArr); + tKine->GetEntry(eventID); + for (int i = 0; i < mcArr->size(); ++i) { + const auto& mcTrack = (*mcArr)[i]; + if (i == trackID) { + Printf("Particle %d: pdg = %d, pT = %f, px = %f, py = %f, pz = %f, vx = %f, vy = %f, vz = %f", i, mcTrack.GetPdgCode(), TMath::Abs(mcTrack.GetStartVertexMomentumX() * mcTrack.GetStartVertexMomentumX() + mcTrack.GetStartVertexMomentumY() * mcTrack.GetStartVertexMomentumY()), mcTrack.GetStartVertexMomentumX(), mcTrack.GetStartVertexMomentumY(), mcTrack.GetStartVertexMomentumZ(), mcTrack.GetStartVertexCoordinatesX(), mcTrack.GetStartVertexCoordinatesY(), mcTrack.GetStartVertexCoordinatesZ()); + } + } + + TFile* fmatch = new TFile("o2match_itstpc.root"); + TTree* matchTPCITS = (TTree*)fmatch->Get("matchTPCITS"); + std::vector* mTracksArrayInp = new std::vector; + matchTPCITS->SetBranchAddress("TPCITS", &mTracksArrayInp); + matchTPCITS->GetEntry(eventID); + + // getting the TPC tracks + TFile* ftracksTPC = new TFile("tpctracks.root"); + TTree* tpcTree = (TTree*)ftracksTPC->Get("events"); + std::vector* mTPCTracksArrayInp = new std::vector; + tpcTree->SetBranchAddress("Tracks", &mTPCTracksArrayInp); + o2::dataformats::MCTruthContainer* mcTPC = new o2::dataformats::MCTruthContainer(); + tpcTree->SetBranchAddress("TPCTracksMCTruth", &mcTPC); + tpcTree->GetEntry(eventID); + + for (int i = 0; i < mTracksArrayInp->size(); i++) { + o2::dataformats::TrackTPCITS trackITSTPC = mTracksArrayInp->at(i); + const o2::dataformats::EvIndex& evIdxTPC = trackITSTPC.getRefTPC(); + const auto& labelsTPC = mcTPC->getLabels(evIdxTPC.getIndex()); + for (int ilabel = 0; ilabel < labelsTPC.size(); ilabel++) { + //Printf("TPC label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsTPC[ilabel].getTrackID(), labelsTPC[ilabel].getEventID(), labelsTPC[ilabel].getSourceID()); + if (labelsTPC[ilabel].getTrackID() == trackID && labelsTPC[ilabel].getEventID() == eventID) + Printf("TPC track found"); + } + } + + return; +} diff --git a/Detectors/TOF/reconstruction/include/TOFReconstruction/Clusterer.h b/Detectors/TOF/reconstruction/include/TOFReconstruction/Clusterer.h index 80a47b2b4a467..06de9d5df5424 100644 --- a/Detectors/TOF/reconstruction/include/TOFReconstruction/Clusterer.h +++ b/Detectors/TOF/reconstruction/include/TOFReconstruction/Clusterer.h @@ -34,7 +34,7 @@ class Clusterer using Digit = o2::tof::Digit; public: - Clusterer(); + Clusterer() = default; ~Clusterer() = default; Clusterer(const Clusterer&) = delete; diff --git a/Detectors/TOF/reconstruction/src/Clusterer.cxx b/Detectors/TOF/reconstruction/src/Clusterer.cxx index 363025373b662..327a5f31711ed 100644 --- a/Detectors/TOF/reconstruction/src/Clusterer.cxx +++ b/Detectors/TOF/reconstruction/src/Clusterer.cxx @@ -19,12 +19,6 @@ using namespace o2::tof; -//__________________________________________________ -Clusterer::Clusterer() -{ - - // empty for now -} //__________________________________________________ void Clusterer::process(DataReader& reader, std::vector& clusters, MCLabelContainer const* digitMCTruth) { @@ -131,7 +125,8 @@ void Clusterer::buildCluster(Cluster& c, MCLabelContainer const* digitMCTruth) } c.setMainContributingChannel(mContributingDigit[0]->getChannel()); - c.setTime(mContributingDigit[0]->getTDC() * Geo::TDCBIN); // time in ps (for now we assume it calibrated) + c.setTime(mContributingDigit[0]->getTDC() * Geo::TDCBIN + double(mContributingDigit[0]->getBC() * 25000.)); // time in ps (for now we assume it calibrated) + c.setTot(mContributingDigit[0]->getTOT() * Geo::TOTBIN * 1E-3); // TOT in ns (for now we assume it calibrated) //setL0L1Latency(); // to be filled (maybe) //setDeltaBC(); // to be filled (maybe) @@ -181,15 +176,27 @@ void Clusterer::buildCluster(Cluster& c, MCLabelContainer const* digitMCTruth) // filling the MC labels of this cluster; the first will be those of the main digit; then the others if (digitMCTruth != nullptr) { int lbl = mClsLabels->getIndexedSize(); // this should correspond to the number of digits also; + printf("lbl = %d\n", lbl); for (int i = 0; i < mNumberOfContributingDigits; i++) { + printf("contributing digit = %d\n", i); int digitLabel = mContributingDigit[i]->getLabel(); + printf("digitLabel = %d\n", digitLabel); gsl::span mcArray = digitMCTruth->getLabels(digitLabel); - for (int j = 0; j < static_cast(mcArray.size()); ++j) { + for (int j = 0; j < static_cast(mcArray.size()); j++) { + printf("checking element %d in the array of labels\n", j); auto label = digitMCTruth->getElement(digitMCTruth->getMCTruthHeader(digitLabel).index + j); + printf("EventID = %d\n", label.getEventID()); mClsLabels->addElement(lbl, label); } } } + // set geometrical variables + int det[5]; + Geo::getVolumeIndices(c.getMainContributingChannel(), det); + float pos[3]; + Geo::getPos(det, pos); + c.setBaseData(c.getMainContributingChannel(), pos[0], pos[1], pos[2], 0, 0, 0); // error on position set to zero + return; } diff --git a/Detectors/TOF/simulation/include/TOFSimulation/Digitizer.h b/Detectors/TOF/simulation/include/TOFSimulation/Digitizer.h index 62dc79679ee3e..9963cc4949583 100644 --- a/Detectors/TOF/simulation/include/TOFSimulation/Digitizer.h +++ b/Detectors/TOF/simulation/include/TOFSimulation/Digitizer.h @@ -22,10 +22,11 @@ namespace o2 { namespace tof { + class Digitizer { public: - Digitizer(Int_t mode = 0) : mMode(mode), mReadoutWindowCurrent(0) { init(); }; + Digitizer(Int_t mode = 0) : mMode(mode) { init(); }; ~Digitizer() = default; void init(); @@ -66,6 +67,9 @@ class Digitizer void setContinuous(bool val) { mContinuous = val; } bool isContinuous() const { return mContinuous; } + std::vector>* getDigitPerTimeFrame() { return &mDigitsPerTimeFrame; } + std::vector>* getMCTruthPerTimeFrame() { return &mMCTruthOutputContainerPerTimeFrame; } + private: // parameters Int_t mMode; @@ -86,7 +90,7 @@ class Digitizer Float_t mEffBoundary3; // info TOF timewindow - Int_t mReadoutWindowCurrent; + Int_t mReadoutWindowCurrent = 0; Double_t mEventTime; Int_t mEventID = 0; Int_t mSrcID = 0; @@ -96,7 +100,10 @@ class Digitizer // digit info //std::vector* mDigits; - static const int MAXWINDOWS = 10; // how many readout windows we can buffer + static const int MAXWINDOWS = 2; // how many readout windows we can buffer + + std::vector> mDigitsPerTimeFrame; + std::vector> mMCTruthOutputContainerPerTimeFrame; int mIcurrentReadoutWindow = 0; o2::dataformats::MCTruthContainer mMCTruthContainer[MAXWINDOWS]; @@ -117,12 +124,14 @@ class Digitizer o2::dataformats::MCTruthContainer mFutureMCTruthContainer; - void fillDigitsInStrip(std::vector* strips, o2::dataformats::MCTruthContainer* mcTruthContainer, double time, int channel, int tdc, int tot, int nbc, UInt_t istrip, Int_t trackID, Int_t eventID, Int_t sourceID); + void fillDigitsInStrip(std::vector* strips, o2::dataformats::MCTruthContainer* mcTruthContainer, int channel, int tdc, int tot, int nbc, UInt_t istrip, Int_t trackID, Int_t eventID, Int_t sourceID); Int_t processHit(const HitType& hit, Double_t event_time); void addDigit(Int_t channel, UInt_t istrip, Float_t time, Float_t x, Float_t z, Float_t charge, Int_t iX, Int_t iZ, Int_t padZfired, Int_t trackID); + void checkIfReuseFutureDigits(); + bool isMergable(Digit digit1, Digit digit2) { if (digit1.getChannel() != digit2.getChannel()) { diff --git a/Detectors/TOF/simulation/include/TOFSimulation/Strip.h b/Detectors/TOF/simulation/include/TOFSimulation/Strip.h index 244991832e500..557854a220b2b 100644 --- a/Detectors/TOF/simulation/include/TOFSimulation/Strip.h +++ b/Detectors/TOF/simulation/include/TOFSimulation/Strip.h @@ -93,7 +93,7 @@ class Strip /// @return Hit at given index (nullptr if index is out of bounds) inline const o2::tof::HitType* getHitAt(Int_t index) const { return mHits.at(index); } - Int_t addDigit(Double_t time, Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t lbl); // returns the MC label + Int_t addDigit(Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t lbl); // returns the MC label void fillOutputContainer(std::vector& digits); diff --git a/Detectors/TOF/simulation/src/Digitizer.cxx b/Detectors/TOF/simulation/src/Digitizer.cxx index e991c59dfc2e5..2f308fb31da86 100644 --- a/Detectors/TOF/simulation/src/Digitizer.cxx +++ b/Detectors/TOF/simulation/src/Digitizer.cxx @@ -49,63 +49,20 @@ void Digitizer::init() void Digitizer::process(const std::vector* hits, std::vector* digits) { // hits array of TOF hits for a given simulated event - //mDigits = digits; + // digits passed from external to be filled, in continuous readout mode we will push it on mDigitsPerTimeFrame vector of vectors of digits - Int_t readoutwindow = Int_t((mEventTime)*Geo::READOUTWINDOW_INV); // to be replaced with uncalibrated time - //while(mContinuous && readoutwindow > mReadoutWindowCurrent){ - if (mContinuous && readoutwindow > mReadoutWindowCurrent) { - digits->clear(); - //fillOutputContainer(*digits); - // waiting for the new framework to store previous readout window + Int_t readoutwindow = Int_t((mEventTime)*Geo::READOUTWINDOW_INV); // to be replaced with "uncalibrated" time - for (Int_t i = mReadoutWindowCurrent; i < readoutwindow; i++) { // temporary loop because current framework doesn't allow to store outputs in digitizer class - fillOutputContainer(*digits); - } - //mReadoutWindowCurrent++; - mReadoutWindowCurrent = readoutwindow; - - // check if digits stored for future match the new readout windows available - int idigit = 0; - for (auto& digit : mFutureDigits) { - int isnext = Int_t(digit.getTimeStamp() * Geo::READOUTWINDOW_INV) - mReadoutWindowCurrent; // to be replaced with uncalibrated time - - if (isnext < 0) // we jump too ahead in future, digit will be not stored - LOG(INFO) << "Digit lost because we jump too ahead in future. Current RO window=" << isnext << "\n"; - - if (isnext < MAXWINDOWS - 1) { // move from digit buffer array to the proper window - if (isnext >= 0) { - std::vector* strips = mStripsCurrent; - o2::dataformats::MCTruthContainer* mcTruthContainer = mMCTruthContainerCurrent; - - if (isnext) { - strips = mStripsNext[isnext - 1]; - mcTruthContainer = mMCTruthContainerNext[isnext - 1]; - } - - int trackID = mFutureItrackID[digit.getLabel()]; - int sourceID = mFutureIsource[digit.getLabel()]; - int eventID = mFutureIevent[digit.getLabel()]; - fillDigitsInStrip(strips, mcTruthContainer, digit.getTimeStamp(), digit.getChannel(), digit.getTDC(), digit.getTOT(), digit.getBC(), digit.getChannel() / Geo::NPADS, trackID, eventID, sourceID); - } + printf("process TOF -> continuous = %i, %i > %i?\n", mContinuous, readoutwindow, mReadoutWindowCurrent); - // remove the element from the buffers - mFutureItrackID.erase(mFutureItrackID.begin() + digit.getLabel()); - mFutureIsource.erase(mFutureIsource.begin() + digit.getLabel()); - mFutureIevent.erase(mFutureIevent.begin() + digit.getLabel()); + if (mContinuous && readoutwindow > mReadoutWindowCurrent) { // if we are moving in future readout windows flush previous ones (only for continuous readout mode) + digits->clear(); - int labelremoved = digit.getLabel(); - // adjust labels - for (auto& digit2 : mFutureDigits) { - if (digit2.getLabel() > labelremoved) - digit2.setLabel(digit2.getLabel() - 1); - } - // remove also digit from buffer - mFutureDigits.erase(mFutureDigits.begin() + idigit); - } else { - idigit++; // increment when moving to the next only if the current is not removed from the buffer - } - } - } + for (; mReadoutWindowCurrent < readoutwindow; mReadoutWindowCurrent++) { + fillOutputContainer(*digits); // fill all windows which are before (not yet stored) of the new current one + checkIfReuseFutureDigits(); + } // close loop readout window + } // close if continuous for (auto& hit : *hits) { //TODO: put readout window counting/selection @@ -294,7 +251,7 @@ void Digitizer::addDigit(Int_t channel, UInt_t istrip, Float_t time, Float_t x, mFutureItrackID.push_back(trackID); // fill temporary digits array - mFutureDigits.emplace_back(time, channel, tdc, tot * Geo::NTOTBIN_PER_NS, nbc, lblCurrent); + mFutureDigits.emplace_back(channel, tdc, tot * Geo::NTOTBIN_PER_NS, nbc, lblCurrent); return; // don't fill if doesn't match any available readout window } @@ -303,6 +260,8 @@ void Digitizer::addDigit(Int_t channel, UInt_t istrip, Float_t time, Float_t x, iscurrent = false; } + // printf("add TOF digit c=%i n=%i\n",iscurrent,isnext); + std::vector* strips; o2::dataformats::MCTruthContainer* mcTruthContainer; @@ -314,17 +273,17 @@ void Digitizer::addDigit(Int_t channel, UInt_t istrip, Float_t time, Float_t x, mcTruthContainer = mMCTruthContainerNext[isnext - 1]; } - fillDigitsInStrip(strips, mcTruthContainer, time, channel, tdc, tot, nbc, istrip, trackID, mEventID, mSrcID); + fillDigitsInStrip(strips, mcTruthContainer, channel, tdc, tot, nbc, istrip, trackID, mEventID, mSrcID); } //______________________________________________________________________ -void Digitizer::fillDigitsInStrip(std::vector* strips, o2::dataformats::MCTruthContainer* mcTruthContainer, double time, int channel, int tdc, int tot, int nbc, UInt_t istrip, Int_t trackID, Int_t eventID, Int_t sourceID) +void Digitizer::fillDigitsInStrip(std::vector* strips, o2::dataformats::MCTruthContainer* mcTruthContainer, int channel, int tdc, int tot, int nbc, UInt_t istrip, Int_t trackID, Int_t eventID, Int_t sourceID) { int lblCurrent; if (mcTruthContainer) { lblCurrent = mcTruthContainer->getIndexedSize(); // this is the size of mHeaderArray; } - Int_t lbl = (*strips)[istrip].addDigit(time, channel, tdc, tot * Geo::NTOTBIN_PER_NS, nbc, lblCurrent); + Int_t lbl = (*strips)[istrip].addDigit(channel, tdc, tot * Geo::NTOTBIN_PER_NS, nbc, lblCurrent); if (mcTruthContainer) { if (lbl == lblCurrent) { // it means that the digit was a new one --> we have to add the info in the MC container @@ -721,12 +680,23 @@ void Digitizer::testFromHits(const char* geo, const char* hits) //______________________________________________________________________ void Digitizer::fillOutputContainer(std::vector& digits) { + if (mContinuous) { + digits.clear(); + mMCTruthOutputContainer->clear(); + } + printf("TOF fill output contatiner\n"); // filling the digit container doing a loop on all strips for (auto& strip : *mStripsCurrent) { strip.fillOutputContainer(digits); } + if (mContinuous) { + if (digits.size()) + printf("%i) # TOF digits = %lu (%p)\n", mIcurrentReadoutWindow, digits.size(), mStripsCurrent); + mDigitsPerTimeFrame.push_back(digits); + } + // if(! digits.size()) return; // copying the transient labels to the output labels (stripping the tdc information) @@ -738,6 +708,8 @@ void Digitizer::fillOutputContainer(std::vector& digits) } } + if (mContinuous) + mMCTruthOutputContainerPerTimeFrame.push_back(*mMCTruthOutputContainer); mMCTruthContainerCurrent->clear(); // switch to next mStrip after flushing current readout window data @@ -745,6 +717,9 @@ void Digitizer::fillOutputContainer(std::vector& digits) if (mIcurrentReadoutWindow >= MAXWINDOWS) mIcurrentReadoutWindow = 0; + mStripsCurrent = &(mStrips[mIcurrentReadoutWindow]); + mMCTruthContainerCurrent = &(mMCTruthContainer[mIcurrentReadoutWindow]); + int k = mIcurrentReadoutWindow + 1; for (Int_t i = 0; i < MAXWINDOWS - 1; i++) { if (k >= MAXWINDOWS) @@ -758,5 +733,64 @@ void Digitizer::fillOutputContainer(std::vector& digits) void Digitizer::flushOutputContainer(std::vector& digits) { // flush all residual buffered data // TO be implemented - fillOutputContainer(digits); + printf("flushOutputContainer\n"); + if (!mContinuous) + fillOutputContainer(digits); + else { + for (Int_t i = 0; i < MAXWINDOWS; i++) { + fillOutputContainer(digits); // fill all windows which are before (not yet stored) of the new current one + checkIfReuseFutureDigits(); + mReadoutWindowCurrent++; + } + + while (mFutureDigits.size()) { + fillOutputContainer(digits); // fill all windows which are before (not yet stored) of the new current one + checkIfReuseFutureDigits(); + mReadoutWindowCurrent++; + } + } +} +//______________________________________________________________________ +void Digitizer::checkIfReuseFutureDigits() +{ + // check if digits stored very far in future match the new readout windows currently available + int idigit = 0; + for (auto& digit : mFutureDigits) { + double timestamp = digit.getBC() * 25 + digit.getTDC() * Geo::TDCBIN * 1E-3; // in ns + int isnext = Int_t(timestamp * Geo::READOUTWINDOW_INV) - (mReadoutWindowCurrent + 1); // to be replaced with uncalibrated time + if (isnext < 0) // we jump too ahead in future, digit will be not stored + LOG(INFO) << "Digit lost because we jump too ahead in future. Current RO window=" << isnext << "\n"; + if (isnext < MAXWINDOWS - 1) { // move from digit buffer array to the proper window + if (isnext >= 0) { + std::vector* strips = mStripsCurrent; + o2::dataformats::MCTruthContainer* mcTruthContainer = mMCTruthContainerCurrent; + + if (isnext) { + strips = mStripsNext[isnext - 1]; + mcTruthContainer = mMCTruthContainerNext[isnext - 1]; + } + + int trackID = mFutureItrackID[digit.getLabel()]; + int sourceID = mFutureIsource[digit.getLabel()]; + int eventID = mFutureIevent[digit.getLabel()]; + fillDigitsInStrip(strips, mcTruthContainer, digit.getChannel(), digit.getTDC(), digit.getTOT(), digit.getBC(), digit.getChannel() / Geo::NPADS, trackID, eventID, sourceID); + } + + // remove the element from the buffers + mFutureItrackID.erase(mFutureItrackID.begin() + digit.getLabel()); + mFutureIsource.erase(mFutureIsource.begin() + digit.getLabel()); + mFutureIevent.erase(mFutureIevent.begin() + digit.getLabel()); + + int labelremoved = digit.getLabel(); + // adjust labels + for (auto& digit2 : mFutureDigits) { + if (digit2.getLabel() > labelremoved) + digit2.setLabel(digit2.getLabel() - 1); + } + // remove also digit from buffer + mFutureDigits.erase(mFutureDigits.begin() + idigit); + } else { + idigit++; // increment when moving to the next only if the current is not removed from the buffer + } + } // close future digit loop } diff --git a/Detectors/TOF/simulation/src/Strip.cxx b/Detectors/TOF/simulation/src/Strip.cxx index 0dccfcacbb4d2..02b9a55e8f8db 100644 --- a/Detectors/TOF/simulation/src/Strip.cxx +++ b/Detectors/TOF/simulation/src/Strip.cxx @@ -56,7 +56,7 @@ Int_t Strip::getStripIndex(const HitType* hit) return channelID / Geo::NPADS; } //_______________________________________________________________________ -Int_t Strip::addDigit(Double_t time, Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t lbl) +Int_t Strip::addDigit(Int_t channel, Int_t tdc, Int_t tot, Int_t bc, Int_t lbl) { // return the MC label. We pass it also as argument, but it can change in @@ -66,9 +66,9 @@ Int_t Strip::addDigit(Double_t time, Int_t channel, Int_t tdc, Int_t tot, Int_t auto dig = findDigit(key); if (dig) { lbl = dig->getLabel(); // getting the label from the already existing digit - dig->merge(time, tdc, tot); // merging to the existing digit + dig->merge(tdc, tot); // merging to the existing digit } else { - auto digIter = mDigits.emplace(std::make_pair(key, Digit(time, channel, tdc, tot, bc, lbl))); + auto digIter = mDigits.emplace(std::make_pair(key, Digit(channel, tdc, tot, bc, lbl))); } return lbl; diff --git a/Generators/src/GeneratorFactory.cxx b/Generators/src/GeneratorFactory.cxx index 758c9de09cf29..b13474a89f343 100644 --- a/Generators/src/GeneratorFactory.cxx +++ b/Generators/src/GeneratorFactory.cxx @@ -176,6 +176,18 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair auto extgen_ptr = (TGenerator**)gROOT->GetGlobal("__extgen")->GetAddress(); tgen->setTGenerator(*extgen_ptr); primGen->AddGenerator(tgen); + } else if (genconfig.compare("toftest") == 0) { // 1 muon per sector and per module + LOG(INFO) << "Init tof test generator -> 1 muon per sector and per module"; + for (int i = 0; i < 18; i++) { + for (int j = 0; j < 5; j++) { + auto boxGen = new FairBoxGenerator(13, 1); /*protons*/ + boxGen->SetEtaRange(-0.8 + 0.32 * j + 0.15, -0.8 + 0.32 * j + 0.17); + boxGen->SetPRange(9, 10); + boxGen->SetPhiRange(10 + 20. * i - 1, 10 + 20. * i + 1); + boxGen->SetDebug(kTRUE); + primGen->AddGenerator(boxGen); + } + } } else { LOG(FATAL) << "Invalid generator"; } diff --git a/Steer/DigitizerWorkflow/src/TOFClusterizerSpec.cxx b/Steer/DigitizerWorkflow/src/TOFClusterizerSpec.cxx index c7dd523b11fe8..cb341f0d4ba93 100644 --- a/Steer/DigitizerWorkflow/src/TOFClusterizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TOFClusterizerSpec.cxx @@ -47,16 +47,20 @@ class TOFDPLClustererTask return; } // get digit data - auto digits = pc.inputs().get*>("tofdigits"); - auto digitlabels = pc.inputs().get*>("tofdigitlabels"); - mReader.setDigitArray(digits.get()); + auto digits = pc.inputs().get>*>("tofdigits"); + auto digitlabels = pc.inputs().get>*>("tofdigitlabels"); mClusterer.setMCTruthContainer(&mClsLabels); // call actual clustering routine mClustersArray.clear(); mClsLabels.clear(); - mClusterer.process(mReader, mClustersArray, digitlabels.get()); + for (int i = 0; i < digits->size(); i++) { + printf("# TOF readout window for clusterization = %i\n", i); + auto digitsRO = digits->at(i); + mReader.setDigitArray(&digitsRO); + mClusterer.process(mReader, mClustersArray, &(digitlabels->at(i))); + } LOG(INFO) << "TOF CLUSTERER : TRANSFORMED " << digits->size() << " DIGITS TO " << mClustersArray.size() << " CLUSTERS"; diff --git a/Steer/DigitizerWorkflow/src/TOFDigitWriterSpec.cxx b/Steer/DigitizerWorkflow/src/TOFDigitWriterSpec.cxx index e1465ec6dacd2..c72f7682d3a9d 100644 --- a/Steer/DigitizerWorkflow/src/TOFDigitWriterSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TOFDigitWriterSpec.cxx @@ -54,7 +54,7 @@ DataProcessorSpec getTOFDigitWriterSpec() auto outputtree = std::make_shared(treename.c_str(), treename.c_str()); // container for incoming digits - auto digits = std::make_shared>(); + auto digits = std::make_shared>>(); // the callback to be set as hook at stop of processing for the framework auto finishWriting = [outputfile, outputtree]() { @@ -76,7 +76,7 @@ DataProcessorSpec getTOFDigitWriterSpec() } // retrieve the digits from the input - auto indata = pc.inputs().get>("tofdigits"); + auto indata = pc.inputs().get>>("tofdigits"); LOG(INFO) << "RECEIVED DIGITS SIZE " << indata.size(); *digits.get() = std::move(indata); @@ -85,8 +85,10 @@ DataProcessorSpec getTOFDigitWriterSpec() br->Fill(); // retrieve labels from the input - auto labeldata = pc.inputs().get*>("tofdigitlabels"); - LOG(INFO) << "TOF GOT " << labeldata->getNElements() << " LABELS "; + auto labeldata = pc.inputs().get>*>("tofdigitlabels"); + for (int i = 0; i < labeldata->size(); i++) { + LOG(INFO) << "TOF GOT " << labeldata->at(i).getNElements() << " LABELS "; + } auto labeldataraw = labeldata.get(); // connect this to a particular branch auto labelbr = getOrMakeBranch(*outputtree.get(), "TOFDigitMCTruth", &labeldataraw); diff --git a/Steer/DigitizerWorkflow/src/TOFDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/TOFDigitizerSpec.cxx index e252c8b6649ed..d563dd97dffd3 100644 --- a/Steer/DigitizerWorkflow/src/TOFDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TOFDigitizerSpec.cxx @@ -17,6 +17,7 @@ #include "TStopwatch.h" #include "Steer/HitProcessingManager.h" // for RunContext #include "TChain.h" +#include "DetectorsBase/GeometryManager.h" #include "TOFSimulation/Digitizer.h" #include "DataFormatsParameters/GRPObject.h" @@ -113,8 +114,8 @@ DataProcessorSpec getTOFDigitizerSpec(int channel) digits->clear(); digitizer->process(&hits, digits.get()); // copy digits into accumulator - std::copy(digits->begin(), digits->end(), std::back_inserter(*digitsAccum.get())); - labelAccum.mergeAtBack(*labels); + //std::copy(digits->begin(), digits->end(), std::back_inserter(*digitsAccum.get())); + //labelAccum.mergeAtBack(*labels); LOG(INFO) << "Have " << digits->size() << " digits "; } } @@ -124,14 +125,23 @@ DataProcessorSpec getTOFDigitizerSpec(int channel) digitizer->flushOutputContainer(*digits.get()); LOG(INFO) << "FLUSHING LEFTOVER STUFF " << digits->size(); // copy digits into accumulator - std::copy(digits->begin(), digits->end(), std::back_inserter(*digitsAccum.get())); - labelAccum.mergeAtBack(*labels); + //std::copy(digits->begin(), digits->end(), std::back_inserter(*digitsAccum.get())); + //labelAccum.mergeAtBack(*labels); + } + + // temporary accumulate vector of vecotors of digits in a single vector + // to be replace once we will be able to write the vector of vectors as different TTree entries + std::vector>* digitsVectOfVect = digitizer->getDigitPerTimeFrame(); + std::vector>* mcLabVecOfVec = digitizer->getMCTruthPerTimeFrame(); + for (Int_t i = 0; i < digitsVectOfVect->size(); i++) { + std::copy(digitsVectOfVect->at(i).begin(), digitsVectOfVect->at(i).end(), std::back_inserter(*digitsAccum.get())); + labelAccum.mergeAtBack(mcLabVecOfVec->at(i)); } LOG(INFO) << "Have " << labelAccum.getNElements() << " TOF labels "; // here we have all digits and we can send them to consumer (aka snapshot it onto output) - pc.outputs().snapshot(Output{ "TOF", "DIGITS", 0, Lifetime::Timeframe }, *digitsAccum.get()); - pc.outputs().snapshot(Output{ "TOF", "DIGITSMCTR", 0, Lifetime::Timeframe }, labelAccum); + pc.outputs().snapshot(Output{ "TOF", "DIGITS", 0, Lifetime::Timeframe }, *digitsVectOfVect); + pc.outputs().snapshot(Output{ "TOF", "DIGITSMCTR", 0, Lifetime::Timeframe }, *mcLabVecOfVec); LOG(INFO) << "TOF: Sending ROMode= " << roMode << " to GRPUpdater"; pc.outputs().snapshot(Output{ "TOF", "ROMode", 0, Lifetime::Timeframe }, roMode); @@ -158,6 +168,11 @@ DataProcessorSpec getTOFDigitizerSpec(int channel) simChains->back()->AddFile(signalfilename.c_str()); } + // make sure that the geometry is loaded (TODO will this be done centrally?) + if (!gGeoManager) { + o2::Base::GeometryManager::loadGeometry(); + } + // init digitizer digitizer->init(); const bool isContinuous = ctx.options().get("pileup"); diff --git a/cmake/O2Dependencies.cmake b/cmake/O2Dependencies.cmake index 947c36ba969db..638b729a35ae9 100644 --- a/cmake/O2Dependencies.cmake +++ b/cmake/O2Dependencies.cmake @@ -2100,12 +2100,14 @@ o2_define_bucket( data_format_reconstruction_bucket data_format_common_bucket data_format_TPC_bucket + data_format_TOF_bucket its_reconstruction_bucket data_format_itsmft_bucket common_field_bucket detectors_base_bucket its_base_bucket tpc_base_bucket + tof_base_bucket data_parameters_bucket common_utils_bucket common_math_bucket @@ -2117,9 +2119,11 @@ o2_define_bucket( DataFormatsITSMFT DetectorsBase DataFormatsTPC + DataFormatsTOF DataFormatsParameters ITSBase TPCBase + TOFBase CommonUtils MathUtils Field @@ -2136,6 +2140,7 @@ o2_define_bucket( ${CMAKE_SOURCE_DIR}/DataFormats/Detectors/ITSMFT/common/include ${CMAKE_SOURCE_DIR}/Detectors/ITSMFT/ITS/base/include ${CMAKE_SOURCE_DIR}/Detectors/TPC/base/include + ${CMAKE_SOURCE_DIR}/Detectors/TOF/base/include ${CMAKE_SOURCE_DIR}/Detectors/Base/include ${CMAKE_SOURCE_DIR}/Common/Utils/include ${CMAKE_SOURCE_DIR}/Common/MathUtils/include diff --git a/macro/checkTOFMatching.C b/macro/checkTOFMatching.C new file mode 100644 index 0000000000000..fba15edf5b890 --- /dev/null +++ b/macro/checkTOFMatching.C @@ -0,0 +1,180 @@ +void checkTOFMatching() +{ + + // macro to check the matching TOF-ITSTPC tracks + + // getting TOF info + TFile* fmatchTOF = new TFile("o2match_tof.root"); + TTree* matchTOF = (TTree*)fmatchTOF->Get("matchTOF"); + std::vector>* TOFMatchInfo; + TOFMatchInfo = new std::vector>; + matchTOF->SetBranchAddress("TOFMatchInfo", &TOFMatchInfo); + + // getting the ITSTPCtracks + TFile* fmatchITSTPC = new TFile("o2match_itstpc.root"); + TTree* matchTPCITS = (TTree*)fmatchITSTPC->Get("matchTPCITS"); + std::vector* mTracksArrayInp = new std::vector; + matchTPCITS->SetBranchAddress("TPCITS", &mTracksArrayInp); + + // getting the TPC tracks + TFile* ftracksTPC = new TFile("tpctracks.root"); + TTree* tpcTree = (TTree*)ftracksTPC->Get("events"); + std::vector* mTPCTracksArrayInp = new std::vector; + tpcTree->SetBranchAddress("TPCTracks", &mTPCTracksArrayInp); + o2::dataformats::MCTruthContainer* mcTPC = new o2::dataformats::MCTruthContainer(); + tpcTree->SetBranchAddress("TPCTracksMCTruth", &mcTPC); + + // getting the ITS tracks + TFile* ftracksITS = new TFile("o2trac_its.root"); + TTree* itsTree = (TTree*)ftracksITS->Get("o2sim"); + std::vector* mITSTracksArrayInp = new std::vector; + itsTree->SetBranchAddress("ITSTrack", &mITSTracksArrayInp); + o2::dataformats::MCTruthContainer* mcITS = new o2::dataformats::MCTruthContainer(); + itsTree->SetBranchAddress("ITSTrackMCTruth", &mcITS); + + // getting the TOF clusters + TFile* fclustersTOF = new TFile("tofclusters.root"); + TTree* tofClTree = (TTree*)fclustersTOF->Get("o2sim"); + std::vector* mTOFClustersArrayInp = new std::vector; + tofClTree->SetBranchAddress("TOFCluster", &mTOFClustersArrayInp); + o2::dataformats::MCTruthContainer* mcTOF = new o2::dataformats::MCTruthContainer(); + tofClTree->SetBranchAddress("TOFClusterMCTruth", &mcTOF); + + tpcTree->GetEntry(0); + tofClTree->GetEntry(0); + + int nMatches = 0; + int nGoodMatches = 0; + int nBadMatches = 0; + + // now looping over the entries in the matching tree + for (int ientry = 0; ientry < matchTOF->GetEntries(); ientry++) { + matchTOF->GetEvent(ientry); + matchTPCITS->GetEntry(ientry); + // now looping over the matched tracks + nMatches += TOFMatchInfo->size(); + for (int imatch = 0; imatch < TOFMatchInfo->size(); imatch++) { + int indexITSTPCtrack = TOFMatchInfo->at(imatch).first; + o2::dataformats::MatchInfoTOF infoTOF = TOFMatchInfo->at(imatch).second; + int tofClIndex = infoTOF.getTOFClIndex(); + float chi2 = infoTOF.getChi2(); + Printf("\nentry in tree %d, matching %d, indexITSTPCtrack = %d, tofClIndex = %d, chi2 = %f", ientry, imatch, indexITSTPCtrack, tofClIndex, chi2); + + // o2::MCCompLabel label = mcTOF->getElement(mcTOF->getMCTruthHeader(tofClIndex).index); + const auto& labelsTOF = mcTOF->getLabels(tofClIndex); + int trackIdTOF; + int eventIdTOF; + int sourceIdTOF; + for (int ilabel = 0; ilabel < labelsTOF.size(); ilabel++) { + Printf("TOF label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsTOF[ilabel].getTrackID(), labelsTOF[ilabel].getEventID(), labelsTOF[ilabel].getSourceID()); + if (ilabel == 0) { + trackIdTOF = labelsTOF[ilabel].getTrackID(); + eventIdTOF = labelsTOF[ilabel].getEventID(); + sourceIdTOF = labelsTOF[ilabel].getSourceID(); + } + } + o2::tof::Cluster tofCluster = mTOFClustersArrayInp->at(tofClIndex); + int nContributingChannels = tofCluster.getNumOfContributingChannels(); + int mainContributingChannel = tofCluster.getMainContributingChannel(); + Printf("The TOF cluster has %d contributing channels, and the main one is %d", nContributingChannels, mainContributingChannel); + int* indices; + o2::tof::Geo::getVolumeIndices(mainContributingChannel, indices); + Printf("Indices of main contributing channel are %d, %d, %d, %d, %d", indices[0], indices[1], indices[2], indices[3], indices[4]); + int* secondaryContributingChannels = new int[nContributingChannels]; + for (int ich = 1; ich < nContributingChannels; ich++) { + bool isUpLeft = tofCluster.isUpLeftContributing(); + bool isUp = tofCluster.isUpContributing(); + bool isUpRight = tofCluster.isUpRightContributing(); + bool isRight = tofCluster.isRightContributing(); + bool isDownRight = tofCluster.isDownRightContributing(); + bool isDown = tofCluster.isDownContributing(); + bool isDownLeft = tofCluster.isDownLeftContributing(); + bool isLeft = tofCluster.isLeftContributing(); + Printf("isUpLeft = %d, isUp = %d, isUpRight = %d, isRight = %d, isDownRight = %d, isDown = %d, isDownLeft = %d, isLeft = %d", isUpLeft, isUp, isUpRight, isRight, isDownRight, isDown, isDownLeft, isLeft); + int* indexCont = new int(); + indexCont[0] = indices[0]; + indexCont[1] = indices[1]; + indexCont[2] = indices[2]; + indexCont[3] = indices[3]; + indexCont[4] = indices[4]; + if (isDown || isDownRight || isDownLeft) { // decrease padZ + indexCont[3]--; + } + if (isUp || isUpRight || isUpLeft) { // decrease padZ + indexCont[3]++; + } + if (isRight || isDownRight || isUpRight) { // decrease padZ + indexCont[4]++; + } + if (isLeft || isDownLeft || isUpLeft) { // decrease padZ + indexCont[4]--; + } + secondaryContributingChannels[ich - 1] = o2::tof::Geo::getIndex(indexCont); + Printf("secondaryContributingChannels[%d] = %d", ich - 1, secondaryContributingChannels[ich - 1]); + } + + o2::dataformats::TrackTPCITS trackITSTPC = mTracksArrayInp->at(indexITSTPCtrack); + const o2::dataformats::EvIndex& evIdxTPC = trackITSTPC.getRefTPC(); + Printf("matched TPCtrack: eventID = %d, indexID = %d", evIdxTPC.getEvent(), evIdxTPC.getIndex()); + const o2::dataformats::EvIndex& evIdxITS = trackITSTPC.getRefITS(); + Printf("matched ITStrack: eventID = %d, indexID = %d", evIdxITS.getEvent(), evIdxITS.getIndex()); + itsTree->GetEntry(evIdxITS.getEvent()); + + // getting the TPC labels + const auto& labelsTPC = mcTPC->getLabels(evIdxTPC.getIndex()); + for (int ilabel = 0; ilabel < labelsTPC.size(); ilabel++) { + Printf("TPC label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsTPC[ilabel].getTrackID(), labelsTPC[ilabel].getEventID(), labelsTPC[ilabel].getSourceID()); + } + + // getting the ITS labels + const auto& labelsITS = mcITS->getLabels(evIdxITS.getIndex()); + for (int ilabel = 0; ilabel < labelsITS.size(); ilabel++) { + Printf("ITS label %d: trackID = %d, eventID = %d, sourceID = %d", ilabel, labelsITS[ilabel].getTrackID(), labelsITS[ilabel].getEventID(), labelsITS[ilabel].getSourceID()); + } + + bool bMatched = kFALSE; + for (int ilabel = 0; ilabel < labelsTOF.size(); ilabel++) { + if ((labelsTPC[0].getTrackID() == labelsTOF[ilabel].getTrackID() && labelsTPC[0].getEventID() == labelsTOF[ilabel].getEventID() && labelsTPC[0].getSourceID() == labelsTOF[ilabel].getSourceID()) || (labelsITS[0].getTrackID() == labelsTOF[ilabel].getTrackID() && labelsITS[0].getEventID() == labelsTOF[ilabel].getEventID() && labelsITS[0].getSourceID() == labelsTOF[ilabel].getSourceID())) { + nGoodMatches++; + bMatched = kTRUE; + break; + } + } + if (!bMatched) + nBadMatches++; + + bool TPCfound = false; + bool ITSfound = false; + for (int i = 0; i < mTracksArrayInp->size(); i++) { + o2::dataformats::TrackTPCITS trackITSTPC = mTracksArrayInp->at(i); + const o2::dataformats::EvIndex& evIdxTPCcheck = trackITSTPC.getRefTPC(); + const o2::dataformats::EvIndex& evIdxITScheck = trackITSTPC.getRefITS(); + itsTree->GetEntry(evIdxITScheck.getEvent()); + const auto& labelsTPCcheck = mcTPC->getLabels(evIdxTPCcheck.getIndex()); + for (int ilabel = 0; ilabel < labelsTPCcheck.size(); ilabel++) { + if (labelsTPCcheck[ilabel].getTrackID() == trackIdTOF && labelsTPCcheck[ilabel].getEventID() == eventIdTOF && labelsTPCcheck[ilabel].getSourceID() == sourceIdTOF) { + Printf("The TPC track that should have been matched to TOF is number %d", i); + TPCfound = true; + } + } + const auto& labelsITScheck = mcITS->getLabels(evIdxITScheck.getIndex()); + for (int ilabel = 0; ilabel < labelsITScheck.size(); ilabel++) { + if (labelsITScheck[ilabel].getTrackID() == trackIdTOF && labelsITScheck[ilabel].getEventID() == eventIdTOF && labelsITScheck[ilabel].getSourceID() == sourceIdTOF) { + Printf("The ITS track that should have been matched to TOF is number %d", i); + ITSfound = true; + } + } + } + if (!TPCfound) + Printf("There is no TPC track found that should have corresponded to this TOF cluster!"); + if (!ITSfound) + Printf("There is no ITS track found that should have corresponded to this TOF cluster!"); + } + } + + Printf("Number of matches = %d", nMatches); + Printf("Number of GOOD matches = %d (%.2f)", nGoodMatches, (float)nGoodMatches / nMatches); + Printf("Number of BAD matches = %d (%.2f)", nBadMatches, (float)nBadMatches / nMatches); + + return; +} diff --git a/macro/run_match_tof.C b/macro/run_match_tof.C new file mode 100644 index 0000000000000..e7adf4799e628 --- /dev/null +++ b/macro/run_match_tof.C @@ -0,0 +1,69 @@ +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include + +#include "Field/MagneticField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" + +#include "GlobalTracking/MatchTOF.h" +#endif + +#define _ALLOW_DEBUG_TREES_ // to allow debug and control tree output + +void run_match_tof(std::string path = "./", std::string outputfile = "o2match_tof.root", + std::string inputTracksTPCITS = "o2match_itstpc.root", + std::string inputTracksTPC = "tpctracks.root", + std::string inputClustersTOF = "tofclusters.root", std::string inputGeom = "O2geometry.root", + std::string inputGRP = "o2sim_grp.root") +{ + + o2::globaltracking::MatchTOF matching; + + if (path.back() != '/') { + path += '/'; + } + + //>>>---------- attach input data --------------->>> + TChain tracks("matchTPCITS"); + tracks.AddFile((path + inputTracksTPCITS).data()); + matching.setInputTreeTracks(&tracks); + + TChain tpcTracks("events"); + tpcTracks.AddFile((path + inputTracksTPC).data()); + matching.setInputTreeTPCTracks(&tpcTracks); + matching.setTPCTrackBranchName("TPCTracks"); + + TChain tofClusters("o2sim"); + tofClusters.AddFile((path + inputClustersTOF).data()); + matching.setInputTreeTOFClusters(&tofClusters); + + //<<<---------- attach input data ---------------<<< + + // create/attach output tree + TFile outFile((path + outputfile).data(), "recreate"); + TTree outTree("matchTOF", "Matched TOF-tracks"); + matching.setOutputTree(&outTree); + +#ifdef _ALLOW_DEBUG_TREES_ + matching.setDebugTreeFileName(path + matching.getDebugTreeFileName()); + matching.setDebugFlag(o2::globaltracking::MatchTOF::MatchTreeAll); +#endif + + //-------- init geometry and field --------// + o2::Base::GeometryManager::loadGeometry(path + inputGeom, "FAIRGeom"); + o2::Base::Propagator::initFieldFromGRP(path + inputGRP); + + matching.init(); + + matching.run(); + + outFile.cd(); + outTree.Write(); + outFile.Close(); +}