Add residuals for external detectors

DataFormats/Detectors/TOF/include/DataFormatsTOF/Cluster.h

@@ -49,14 +49,19 @@ class Cluster : public o2::BaseCluster<float>
          kDownRight = 4, // 2^4, 5th bit
          kDown = 5,      // 2^5, 6th bit
          kDownLeft = 6,  // 2^6, 7th bit
-         kLeft = 7 };    // 2^7, 8th bit
+         kLeft = 7,      // 2^7, 8th bit
+         //
+         FrameBit = 6 }; // this bit set means that the cluster is in the nominal (alpha=20*sector+10 deg.) sector frame rather than aligned
 
   Cluster() = default;
 
   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, float geanttime = 0.0, double t0 = 0.0);
 
   ~Cluster() = default;
 
+  bool isInNominalSector() const { return isBitSet(FrameBit); }
+  void setInNominalSector() { setBit(FrameBit); }
+
   std::int8_t getSector() const { return getCount(); }
   void setSector(std::int8_t value) { setCount(value); }
 

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/include/TPCInterpolationWorkflow/TPCInterpolationSpec.h

@@ -38,7 +38,8 @@ namespace tpc
 class TPCInterpolationDPL : public Task
 {
  public:
-  TPCInterpolationDPL(std::shared_ptr<o2::globaltracking::DataRequest> dr, o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap, std::shared_ptr<o2::base::GRPGeomRequest> gr, bool useMC, bool processITSTPConly, bool sendTrackData, bool debugOutput) : mDataRequest(dr), mSources(src), mSourcesMap(srcMap), mGGCCDBRequest(gr), mUseMC(useMC), mProcessITSTPConly(processITSTPConly), mSendTrackData(sendTrackData), mDebugOutput(debugOutput) {}
+  TPCInterpolationDPL(std::shared_ptr<o2::globaltracking::DataRequest> dr, o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap, std::shared_ptr<o2::base::GRPGeomRequest> gr, bool useMC,
+                      bool processITSTPConly, bool sendTrackData, bool debugOutput, bool extDetResid) : mDataRequest(dr), mSources(src), mSourcesMap(srcMap), mGGCCDBRequest(gr), mUseMC(useMC), mProcessITSTPConly(processITSTPConly), mSendTrackData(sendTrackData), mDebugOutput(debugOutput), mExtDetResid(extDetResid) {}
   ~TPCInterpolationDPL() override = default;
   void init(InitContext& ic) final;
   void run(ProcessingContext& pc) final;
@@ -58,14 +59,16 @@ class TPCInterpolationDPL : public Task
   bool mProcessITSTPConly{false}; ///< should also tracks without outer point (ITS-TPC only) be processed?
   bool mProcessSeeds{false};      ///< process not only most complete track, but also its shorter parts
   bool mDebugOutput{false};       ///< add more information to the output (track points of ITS, TRD and TOF)
+  bool mExtDetResid{true};        ///< produce unbinned residuals for external detectors
   bool mSendTrackData{false};     ///< if true, not only the clusters but also corresponding track data will be sent
   uint32_t mSlotLength{600u};     ///< the length of one calibration slot required to calculate max number of tracks per TF
   int mMatCorr{2};                ///< the material correction to be used for track interpolation
   TStopwatch mTimer;
 };
 
 /// create a processor spec
-framework::DataProcessorSpec getTPCInterpolationSpec(o2::dataformats::GlobalTrackID::mask_t srcCls, o2::dataformats::GlobalTrackID::mask_t srcVtx, o2::dataformats::GlobalTrackID::mask_t srcTrk, o2::dataformats::GlobalTrackID::mask_t srcTrkMap, bool useMC, bool processITSTPConly, bool sendTrackData, bool debugOutput);
+framework::DataProcessorSpec getTPCInterpolationSpec(o2::dataformats::GlobalTrackID::mask_t srcCls, o2::dataformats::GlobalTrackID::mask_t srcVtx, o2::dataformats::GlobalTrackID::mask_t srcTrk,
+                                                     o2::dataformats::GlobalTrackID::mask_t srcTrkMap, bool useMC, bool processITSTPConly, bool sendTrackData, bool debugOutput, bool extDetResid);
 
 } // namespace tpc
 } // namespace o2

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx

@@ -105,8 +105,9 @@ void TPCInterpolationDPL::updateTimeDependentParams(ProcessingContext& pc)
   }
   if (mDebugOutput) {
     mInterpolation.setDumpTrackPoints();
-    mInterpolation.setITSClusterDictionary(mITSDict);
   }
+  mInterpolation.setExtDetResid(mExtDetResid);
+  mInterpolation.setITSClusterDictionary(mITSDict);
 }
 
 void TPCInterpolationDPL::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
@@ -158,7 +159,7 @@ void TPCInterpolationDPL::endOfStream(EndOfStreamContext& ec)
        mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1);
 }
 
-DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mask_t srcVtx, GTrackID::mask_t srcTrk, GTrackID::mask_t srcTrkMap, bool useMC, bool processITSTPConly, bool sendTrackData, bool debugOutput)
+DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mask_t srcVtx, GTrackID::mask_t srcTrk, GTrackID::mask_t srcTrkMap, bool useMC, bool processITSTPConly, bool sendTrackData, bool debugOutput, bool extDetResid)
 {
   auto dataRequest = std::make_shared<DataRequest>();
   std::vector<OutputSpec> outputs;
@@ -199,7 +200,7 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas
     "tpc-track-interpolation",
     dataRequest->inputs,
     outputs,
-    AlgorithmSpec{adaptFromTask<TPCInterpolationDPL>(dataRequest, srcTrk, srcTrkMap, ggRequest, useMC, processITSTPConly, sendTrackData, debugOutput)},
+    AlgorithmSpec{adaptFromTask<TPCInterpolationDPL>(dataRequest, srcTrk, srcTrkMap, ggRequest, useMC, processITSTPConly, sendTrackData, debugOutput, extDetResid)},
     Options{
       {"matCorrType", VariantType::Int, 2, {"material correction type (definition in Propagator.h)"}},
       {"sec-per-slot", VariantType::UInt32, 600u, {"number of seconds per calibration time slot (put 0 for infinite slot length)"}},

Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/tpc-interpolation-workflow.cxx

@@ -42,6 +42,7 @@ void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
     {"tracking-sources-map-extraction", VariantType::String, std::string{GID::ALL}, {"can be subset of \"tracking-sources\""}},
     {"send-track-data", VariantType::Bool, false, {"Send also the track information to the aggregator"}},
     {"debug-output", VariantType::Bool, false, {"Dump extended tracking information for debugging"}},
+    {"skip-ext-det-residuals", VariantType::Bool, false, {"Do not produce residuals for external detectors"}},
     {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}};
   o2::raw::HBFUtilsInitializer::addConfigOption(options);
   std::swap(workflowOptions, options);
@@ -104,8 +105,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext)
   useMC = false; // force disabling MC as long as it is not implemented
   auto sendTrackData = configcontext.options().get<bool>("send-track-data");
   auto debugOutput = configcontext.options().get<bool>("debug-output");
+  auto extDetResid = !configcontext.options().get<bool>("skip-ext-det-residuals");
 
-  specs.emplace_back(o2::tpc::getTPCInterpolationSpec(srcClusters, srcVtx, srcTracks, srcTracksMap, useMC, processITSTPConly, sendTrackData, debugOutput));
+  specs.emplace_back(o2::tpc::getTPCInterpolationSpec(srcClusters, srcVtx, srcTracks, srcTracksMap, useMC, processITSTPConly, sendTrackData, debugOutput, extDetResid));
   if (!configcontext.options().get<bool>("disable-root-output")) {
     specs.emplace_back(o2::tpc::getTPCResidualWriterSpec(sendTrackData, debugOutput));
   }

Detectors/TOF/base/include/TOFBase/Geo.h

@@ -66,6 +66,8 @@ class Geo
   static void antiRotateToSector(Float_t* xyz, Int_t isector);
   static void antiRotateToStrip(Float_t* xyz, Int_t iplate, Int_t istrip, Int_t isector);
 
+  static void alignedToNominalSector(Float_t* xyz, Int_t isector);
+
   static void antiRotate(Float_t* xyz, Double_t rotationAngles[6]);
   static void getDetID(Float_t* pos, Int_t* det);
   static Int_t getIndex(const Int_t* detId);               // Get channel index from det Id (for calibration mainly)

Detectors/TOF/base/src/Geo.cxx

@@ -988,6 +988,7 @@ void Geo::translate(Float_t* xyz, Float_t translationVector[3])
 
   return;
 }
+
 void Geo::translate(Float_t& x, Float_t& y, Float_t& z, Float_t translationVector[3])
 {
   //
@@ -1045,6 +1046,18 @@ void Geo::rotateToSector(Float_t* xyz, Int_t isector)
   return;
 }
 
+void Geo::alignedToNominalSector(Float_t* xyz, Int_t isector)
+{
+  // rotate from the aligned sector frame coordinates to nominal ones (i.e. alpha=20*sector+10 deg.)
+  constexpr float CS[18] = {.848077e-01, 8.660254e-01, 6.427876e-01, 3.420202e-01, -4.371139e-08, -3.420201e-01, -6.427876e-01, -8.660254e-01, -9.848077e-01, -9.848077e-01, -8.660254e-01, -6.427875e-01, -3.420201e-01, 1.192488e-08, 3.420201e-01, 6.427875e-01, 8.660253e-01, 9.848078e-01};
+  constexpr float SN[18] = {1.736482e-01, 5.000000e-01, 7.660444e-01, 9.396926e-01, 1.000000e+00, 9.396926e-01, 7.660444e-01, 5.000001e-01, 1.736483e-01, -1.736482e-01, -5.000000e-01, -7.660446e-01, -9.396927e-01, -1.000000e+00, -9.396926e-01, -7.660445e-01, -5.000002e-01, -1.736480e-01};
+  Float_t xyzDummy[3] = {xyz[1], xyz[2], xyz[0]};      // go to twisted coordinates...
+  o2::tof::Geo::antiRotateToSector(xyzDummy, isector); // lab coordinates
+  xyz[0] = xyzDummy[0] * CS[isector] + xyzDummy[1] * SN[isector];
+  xyz[1] = -xyzDummy[0] * SN[isector] + xyzDummy[1] * CS[isector];
+  xyz[2] = xyzDummy[2];
+}
+
 void Geo::antiRotateToStrip(Float_t* xyz, Int_t iplate, Int_t istrip, Int_t isector)
 {
   Float_t xyzDummy[3] = {0., 0., 0.};

Detectors/TPC/calibration/SpacePoints/CMakeLists.txt

@@ -19,6 +19,7 @@ o2_add_library(SpacePoints
                                      O2::CommonUtils
                                      O2::TPCBase
                                      O2::TRDBase
+                                     O2::TOFBase
                                      O2::TPCReconstruction
                                      O2::TPCFastTransformation
                                      O2::ITStracking

Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h

@@ -67,31 +67,49 @@ struct TPCClusterResiduals {
 /// (this is the data type which will be sent from the EPNs to the aggregator)
 struct UnbinnedResid {
   UnbinnedResid() = default;
-  UnbinnedResid(float dyIn, float dzIn, float tgSlpIn, float yIn, float zIn, unsigned char rowIn, unsigned char secIn) : dy(static_cast<short>(dyIn * 0x7fff / param::MaxResid)),
-                                                                                                                         dz(static_cast<short>(dzIn * 0x7fff / param::MaxResid)),
-                                                                                                                         tgSlp(static_cast<short>(tgSlpIn * 0x7fff / param::MaxTgSlp)),
-                                                                                                                         y(static_cast<short>(yIn * 0x7fff / param::MaxY)),
-                                                                                                                         z(static_cast<short>(zIn * 0x7fff / param::MaxZ)),
-                                                                                                                         row(rowIn),
-                                                                                                                         sec(secIn) {}
-  short dy;          ///< residual in y
-  short dz;          ///< residual in z
-  short tgSlp;       ///< tan of the phi angle between padrow and track
-  short y;           ///< y position of the track, needed for binning
-  short z;           ///< z position of the track, needed for binning
-  unsigned char row; ///< TPC pad row
-  unsigned char sec; ///< TPC sector (0..35)
-  ClassDefNV(UnbinnedResid, 1);
+  UnbinnedResid(float dyIn, float dzIn, float tgSlpIn, float yIn, float zIn, unsigned char rowIn, unsigned char secIn, short chanIn = -1) : dy(static_cast<short>(dyIn * 0x7fff / param::MaxResid)),
+                                                                                                                                            dz(static_cast<short>(dzIn * 0x7fff / param::MaxResid)),
+                                                                                                                                            tgSlp(static_cast<short>(tgSlpIn * 0x7fff / param::MaxTgSlp)),
+                                                                                                                                            y(static_cast<short>(yIn * 0x7fff / param::MaxY)),
+                                                                                                                                            z(static_cast<short>(zIn * 0x7fff / param::MaxZ)),
+                                                                                                                                            row(rowIn),
+                                                                                                                                            sec(secIn),
+                                                                                                                                            channel(chanIn) {}
+  short dy{0};          ///< residual in y
+  short dz{0};          ///< residual in z
+  short tgSlp{0};       ///< tan of the phi angle between padrow and track
+  short y{0};           ///< y position of the track, needed for binning
+  short z{0};           ///< z position of the track, needed for binning
+  unsigned char row{0}; ///< TPC pad row
+  unsigned char sec{0}; ///< TPC sector (0..35)
+  short channel{-1};    ///< extra channel info (ITS chip ID, TRD chamber, TOF main pad within the sector)
+
+  bool isTPC() const { return row < constants::MAXGLOBALPADROW; }
+  bool isTRD() const { return row >= 160 && row < 166; }
+  bool isTOF() const { return row == 170; }
+  bool isITS() const { return row >= 180; }
+  int getDetID() const { return isTPC() ? 1 : (isITS() ? 0 : (isTRD() ? 2 : (isTOF() ? 3 : -1))); }
+  int getITSLayer() const { return row - 180; }
+  int getTRDLayer() const { return row - 170; }
+  float getAlpha() const;
+  float getX() const;
+
+  static void init(long timestamp = -1);
+  static void checkInitDone();
+  static bool gInitDone;
+
+  ClassDefNV(UnbinnedResid, 2);
 };
 
 /// Structure for the information required to associate each residual with a given track type (ITS-TPC-TRD-TOF, etc)
 struct TrackDataCompact {
   TrackDataCompact() = default;
-  TrackDataCompact(uint32_t idx, uint8_t nRes, uint8_t source) : idxFirstResidual(idx), nResiduals(nRes), sourceId(source) {}
+  TrackDataCompact(uint32_t idx, uint8_t nRes, uint8_t source, uint8_t nextraRes = 0) : idxFirstResidual(idx), nResiduals(nRes), sourceId(source), nExtDetResid(nextraRes) {}
   uint32_t idxFirstResidual; ///< the index of the first residual from this track
-  uint8_t nResiduals;        ///< total number of residuals associated to this track
+  uint8_t nResiduals;        ///< total number of TPC residuals associated to this track
+  uint8_t nExtDetResid = 0;  ///< number of external detectors (wrt TPC) residuals stored, on top of clIdx.getEntries
   uint8_t sourceId;          ///< source ID obtained from the global track ID
-  ClassDefNV(TrackDataCompact, 1);
+  ClassDefNV(TrackDataCompact, 2);
 };
 
 // TODO add to UnbinnedResid::sec flag if cluster was used or not
@@ -110,7 +128,8 @@ struct TrackDataExtended {
   std::vector<o2::trd::CalibratedTracklet> clsTRD{}; ///< the TRD space points (if available)
   o2::tof::Cluster clsTOF{};                         ///< the TOF cluster (if available)
   o2::dataformats::RangeReference<> clIdx{};         ///< index of first cluster residual and total number of cluster residuals of this track
-  ClassDefNV(TrackDataExtended, 2);
+  uint8_t nExtDetResid = 0;                          ///< number of external detectors (to TPC) residuals stored, on top of clIdx.getEntries
+  ClassDefNV(TrackDataExtended, 3);
 };
 
 /// Structure filled for each track with track quality information and a vector with TPCClusterResiduals
@@ -121,12 +140,14 @@ struct TrackData {
   float chi2TPC{};                           ///< chi2 of TPC track
   float chi2ITS{};                           ///< chi2 of ITS track
   float chi2TRD{};                           ///< chi2 of TRD track
+
   unsigned short nClsTPC{};                  ///< number of attached TPC clusters
   unsigned short nClsITS{};                  ///< number of attached ITS clusters
   unsigned short nTrkltsTRD{};               ///< number of attached TRD tracklets
   unsigned short clAvailTOF{};               ///< whether or not track seed has a matched TOF cluster
-  o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of cluster residuals of this track
-  ClassDefNV(TrackData, 6);
+  uint8_t nExtDetResid = 0;                  ///< number of external detectors (to TPC) residuals stored, on top of clIdx.getEntries
+  o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of TPC cluster residuals of this track
+  ClassDefNV(TrackData, 7);
 };
 
 /// \class TrackInterpolation
@@ -265,6 +286,10 @@ class TrackInterpolation
   /// Set the centre of mass energy required for pT downsampling Tsalis function
   void setSqrtS(float s) { mSqrtS = s; }
 
+  void setExtDetResid(bool v) { mExtDetResid = v; }
+
+  int processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLayer, o2::track::TrackParCov& trkWork, std::array<float, 2>* trkltTRDYZ = nullptr, std::array<float, 3>* trkltTRDCov = nullptr);
+
   // --------------------------------- output ---------------------------------------------
   std::vector<UnbinnedResid>& getClusterResiduals() { return mClRes; }
   std::vector<TrackDataCompact>& getTrackDataCompact() { return mTrackDataCompact; }
@@ -285,6 +310,7 @@ class TrackInterpolation
   int mMaxTracksPerTF{-1};                                      ///< max number of tracks to be processed per TF (-1 means there is no limit)
   int mAddTracksForMapPerTF{0};                                 ///< in case residuals from different track types are used for vDrift calibration and map creation this defines the statistics for the latter
   bool mDumpTrackPoints{false};                                 ///< dump also track points in ITS, TRD and TOF
+  bool mExtDetResid{true};                                      ///< produce unbinned residuals for external detectors
   bool mProcessSeeds{false};                                    ///< in case for global tracks also their shorter parts are processed separately
   bool mProcessITSTPConly{false};                               ///< flag, whether or not to extrapolate ITS-only through TPC
   o2::dataformats::GlobalTrackID::mask_t mSourcesConfigured;    ///< the track sources taken into account for extra-/interpolation
@@ -297,6 +323,7 @@ class TrackInterpolation
   std::vector<o2::globaltracking::RecoContainer::GlobalIDSet> mGIDtables{}; ///< GIDs of contributors from single detectors for each seed
   std::vector<float> mTrackTimes{};                                         ///< time estimates for all input tracks in micro seconds
   std::vector<o2::track::TrackParCov> mSeeds{};                             ///< seeding track parameters (ITS tracks)
+  std::vector<int> mParentID{};                                             ///< entry of more global parent track for skimmed seeds (-1: no parent)
   std::map<int, int> mTrackTypes;                                           ///< mapping of track source to array index in mTrackIndices
   std::array<std::vector<uint32_t>, 4> mTrackIndices;                       ///< keep GIDs of input tracks separately for each track type
   gsl::span<const TPCClRefElem> mTPCTracksClusIdx;                          ///< input TPC cluster indices from span