Extend TOF Event Time maker for analysis

Detectors/AOD/src/AODProducerWorkflowSpec.cxx

@@ -327,8 +327,11 @@ void AODProducerWorkflowDPL::fillTrackTablesPerCollision(int collisionID,
             }
             const auto& tofCl = tofClus[contributorsGID[GIndex::Source::TOF]];
             // correct the time of the track
-            extraInfoHolder.trackTime = (tofCl.getTime() - tofInt.getTOF(trackPar.getPID())) * 1e-3; // tof time in \mus, FIXME: account for time of flight to R TOF
-            extraInfoHolder.trackTimeRes = 200e-3;                                                   // FIXME: calculate actual resolution (if possible?)
+            const float massZ = o2::track::PID::getMass2Z(trackPar.getPID());
+            const float energy = sqrt((massZ * massZ) + (extraInfoHolder.tofExpMom * extraInfoHolder.tofExpMom));
+            const float exp = extraInfoHolder.length * energy / (cSpeed * extraInfoHolder.tofExpMom);
+            extraInfoHolder.trackTime = (tofCl.getTime() - exp) * 1e-3; // tof time in \mus, FIXME: account for time of flight to R TOF
+            extraInfoHolder.trackTimeRes = 200e-3;                      // FIXME: calculate actual resolution (if possible?)
           }
           if (src == GIndex::Source::TPCTRD || src == GIndex::Source::ITSTPCTRD) {
             const auto& trdOrig = data.getTrack<o2::trd::TrackTRD>(src, contributorsGID[src].getIndex());

Detectors/GlobalTrackingWorkflow/src/TOFEventTimeChecker.cxx

@@ -64,7 +64,7 @@ using GID = o2::dataformats::GlobalTrackID;
 
 struct MyTrack : o2::tof::eventTimeTrackTest {
   double tofSignalDouble() const { return mSignalDouble; }
-  float tofExpTimeDe() const { return mExpDe; }
+  float tofExpSignalDe() const { return mExpDe; }
   double mSignalDouble = 0.0;
   float mEta = 0.0;
   float mPhi = 0.0;
@@ -216,10 +216,10 @@ void TOFEventTimeChecker::processEvent(std::vector<MyTrack>& tracks)
     mChi2 = track.mChi2;
     mL = track.mLength;
     mTof = track.tofSignal();
-    mExpDe = track.tofExpTimeDe();
-    mExpPi = track.tofExpTimePi();
-    mExpKa = track.tofExpTimeKa();
-    mExpPr = track.tofExpTimePr();
+    mExpDe = track.tofExpSignalDe();
+    mExpPi = track.tofExpSignalPi();
+    mExpKa = track.tofExpSignalKa();
+    mExpPr = track.tofExpSignalPr();
     mIsProb = track.mIsProb;
     mTrktime = track.mTrktime;
     mTrktimeRes = track.mTrktimeRes;

Detectors/TOF/prototyping/macroEvTime.C

@@ -45,10 +45,10 @@ void macroEvTime(bool removebias = true)
       nt++;
 
       float res = sqrt(100 * 100 + erret * erret);
-      hkaon->Fill((track.tofSignal() - track.tofExpTimeKa() - et) / res);
+      hkaon->Fill((track.tofSignal() - track.tofExpSignalKa() - et) / res);
 
       if (track.masshypo() == 1) {
-        hkaonT->Fill((track.tofSignal() - track.tofExpTimeKa() - et) / res);
+        hkaonT->Fill((track.tofSignal() - track.tofExpSignalKa() - et) / res);
       }
     }
   }

Detectors/TOF/reconstruction/include/TOFReconstruction/EventTimeMaker.h

@@ -9,15 +9,20 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
-/// \file EventTimeMaker.h
-/// \brief Definition of the TOF event time maker
+/// \file   EventTimeMaker.h
+/// \author Francesca Ercolessi francesca.ercolessi@cern.ch
+/// \author Francesco Noferini francesco.noferini@cern.ch
+/// \author Nicolò Jacazio nicolo.jacazio@cern.ch
+/// \brief  Definition of the TOF event time maker
 
 #ifndef ALICEO2_TOF_EVENTTIMEMAKER_H
 #define ALICEO2_TOF_EVENTTIMEMAKER_H
 
 #include "TRandom.h"
 #include "TMath.h"
 #include "TOFBase/Utils.h"
+#include "ReconstructionDataFormats/PID.h"
+#include "Framework/Logger.h"
 
 namespace o2
 {
@@ -27,12 +32,17 @@ namespace tof
 
 struct eventTimeContainer {
   eventTimeContainer(const float& e, const float& err) : eventTime{e}, eventTimeError{err} {};
-  float eventTime = 0.f;
-  float eventTimeError = 0.f;
-
-  float sumweights = 0.f;       // sum of weights of all track contributors
-  std::vector<float> weights;   // weights (1/sigma^2) associated to a track in event time computation, 0 if track not used
-  std::vector<float> tracktime; // eventtime provided by a single track
+  float eventTime = 0.f;                      /// Value of the event time
+  float eventTimeError = 0.f;                 /// Uncertainty on the computed event time
+  unsigned short eventTimeMultiplicity = 0.f; /// Track multiplicity used to compute the event time
+
+  float sumweights = 0.f;       /// sum of weights of all track contributors
+  std::vector<float> weights;   /// weights (1/sigma^2) associated to a track in event time computation, 0 if track not used
+  std::vector<float> tracktime; /// eventtime provided by a single track
+  void print()
+  {
+    LOG(info) << "eventTimeContainer " << eventTime << " +- " << eventTimeError << " sum of weights " << sumweights << " tracks used " << eventTimeMultiplicity;
+  }
 };
 
 struct eventTimeTrack {
@@ -45,9 +55,9 @@ struct eventTimeTrack {
     }
   }
   float tofSignal() const { return mSignal; }
-  float tofExpTimePi() const { return expTimes[0]; }
-  float tofExpTimeKa() const { return expTimes[1]; }
-  float tofExpTimePr() const { return expTimes[2]; }
+  float tofExpSignalPi() const { return expTimes[0]; }
+  float tofExpSignalKa() const { return expTimes[1]; }
+  float tofExpSignalPr() const { return expTimes[2]; }
   float tofExpSigmaPi() const { return expSigma[0]; }
   float tofExpSigmaKa() const { return expSigma[1]; }
   float tofExpSigmaPr() const { return expSigma[2]; }
@@ -74,14 +84,17 @@ void generateEvTimeTracks(std::vector<eventTimeTrackTest>& tracks, int ntracks,
 template <typename trackType>
 bool filterDummy(const trackType& tr)
 {
-  return (tr.tofChi2() >= 0 && tr.mP < 2.0);
+  return (tr.tofChi2() >= 0 && tr.p() < 2.0);
 } // accept all
 
 void computeEvTime(const std::vector<eventTimeTrack>& tracks, const std::vector<int>& trkIndex, eventTimeContainer& evtime);
 int getStartTimeInSet(const std::vector<eventTimeTrack>& tracks, std::vector<int>& trackInSet, unsigned long& bestComb);
 
-template <typename trackTypeContainer, typename trackType, bool (*trackFilter)(const trackType&)>
-eventTimeContainer evTimeMaker(const trackTypeContainer& tracks, float diamond = Utils::mEventTimeSpread * 0.029979246 /* spread of primary verdex in cm */)
+template <typename trackTypeContainer,
+          typename trackType,
+          bool (*trackFilter)(const trackType&)>
+eventTimeContainer evTimeMaker(const trackTypeContainer& tracks,
+                               float diamond = Utils::mEventTimeSpread * 0.029979246 /* spread of primary verdex in cm */)
 {
   static std::vector<eventTimeTrack> trkWork;
   trkWork.clear();
@@ -100,11 +113,14 @@ eventTimeContainer evTimeMaker(const trackTypeContainer& tracks, float diamond =
   result.eventTimeError = sigmaFill;
   result.sumweights = 0.;
 
-  // Qui facciamo un pool di tracce buone per calcolare il T0
-  for (auto track : tracks) {
-    if (trackFilter(track)) {
-      expt[0] = track.tofExpTimePi(), expt[1] = track.tofExpTimeKa(), expt[2] = track.tofExpTimePr();
-      expsigma[0] = track.tofExpSigmaPi(), expsigma[1] = track.tofExpSigmaKa(), expsigma[2] = track.tofExpSigmaPr();
+  for (auto track : tracks) { // Loop on tracks
+    if (trackFilter(track)) { // Select tracks good for T0 computation
+      expt[0] = track.tofExpSignalPi();
+      expt[1] = track.tofExpSignalKa();
+      expt[2] = track.tofExpSignalPr();
+      expsigma[0] = track.tofExpSigmaPi();
+      expsigma[1] = track.tofExpSigmaKa();
+      expsigma[2] = track.tofExpSigmaPr();
       trkWork.emplace_back(track.tofSignal(), expt, expsigma);
       trkIndex.push_back(result.weights.size());
     }
@@ -115,6 +131,53 @@ eventTimeContainer evTimeMaker(const trackTypeContainer& tracks, float diamond =
   return result;
 }
 
+template <typename trackTypeContainer,
+          typename trackType,
+          bool (*trackFilter)(const trackType&),
+          template <typename T, o2::track::PID::ID> typename response,
+          typename responseParametersType>
+eventTimeContainer evTimeMakerFromParam(const trackTypeContainer& tracks,
+                                        const responseParametersType& responseParameters,
+                                        float diamond = Utils::mEventTimeSpread * 0.029979246 /* spread of primary verdex in cm */)
+{
+  static std::vector<eventTimeTrack> trkWork;
+  trkWork.clear();
+  static std::vector<int> trkIndex; // indexes of working tracks in the track original array
+  trkIndex.clear();
+
+  constexpr auto responsePi = response<trackType, o2::track::PID::Pion>();
+  constexpr auto responseKa = response<trackType, o2::track::PID::Kaon>();
+  constexpr auto responsePr = response<trackType, o2::track::PID::Proton>();
+
+  static float expt[3], expsigma[3];
+
+  static eventTimeContainer result = {0, 0};
+
+  // reset info
+  float sigmaFill = diamond * 33.356409; // move from diamond (cm) to spread on event time (ps)
+  result.weights.clear();
+  result.tracktime.clear();
+  result.eventTime = 0.;
+  result.eventTimeError = sigmaFill;
+  result.sumweights = 0.;
+
+  for (auto track : tracks) { // Loop on tracks
+    if (trackFilter(track)) { // Select tracks good for T0 computation
+      expt[0] = responsePi.GetExpectedSignal(track);
+      expt[1] = responseKa.GetExpectedSignal(track);
+      expt[2] = responsePr.GetExpectedSignal(track);
+      expsigma[0] = responsePi.GetExpectedSigmaTracking(responseParameters, track);
+      expsigma[1] = responseKa.GetExpectedSigmaTracking(responseParameters, track);
+      expsigma[2] = responsePr.GetExpectedSigmaTracking(responseParameters, track);
+      trkWork.emplace_back(track.tofSignal(), expt, expsigma);
+      trkIndex.push_back(result.weights.size());
+    }
+    result.weights.push_back(0.);
+    result.tracktime.push_back(0.);
+  }
+  computeEvTime(trkWork, trkIndex, result);
+  return result;
+}
 } // namespace tof
 } // namespace o2
 

Detectors/TOF/reconstruction/src/EventTimeMaker.cxx

@@ -9,9 +9,16 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
-/// \file EventTimeMaker.cxx
-/// \brief Implementation of the TOF event time maker
-
+///
+/// \file   EventTimeMaker.cxx
+/// \author Francesca Ercolessi francesca.ercolessi@cern.ch
+/// \author Francesco Noferini francesco.noferini@cern.ch
+/// \author Nicolò Jacazio nicolo.jacazio@cern.ch
+/// \brief  Implementation of the TOF event time maker
+///
+
+#include "TRandom.h"
+#include "TMath.h"
 #include "TOFReconstruction/EventTimeMaker.h"
 
 namespace o2
@@ -27,11 +34,14 @@ constexpr unsigned long combinatorial[MAXNTRACKINSET + 1] = {1, 3, 9, 27, 81, 24
 
 void computeEvTime(const std::vector<eventTimeTrack>& tracks, const std::vector<int>& trkIndex, eventTimeContainer& evtime)
 {
-  const int maxNumberOfSets = 200;
+  static constexpr int maxNumberOfSets = 200;
+  static constexpr float weightLimit = 1E-6; // Limit in the weights
 
-  int ntracks = tracks.size();
+  const int ntracks = tracks.size();
+  LOG(debug) << "For the collision time using " << ntracks;
 
   if (ntracks < 2) { // at least 2 tracks required
+    LOG(debug) << "Skipping event because at least 2 tracks are required";
     return;
   }
 
@@ -40,7 +50,7 @@ void computeEvTime(const std::vector<eventTimeTrack>& tracks, const std::vector<
   int nmaxtracksinset = ntracks > 22 ? 6 : MAXNTRACKINSET; // max number of tracks in a set for event time computation
   int ntracksinset = std::min(ntracks, nmaxtracksinset);
 
-  Int_t nset = ((ntracks - 1) / ntracksinset) + 1;
+  int nset = ((ntracks - 1) / ntracksinset) + 1;
   int ntrackUsed = ntracks;
 
   if (nset > maxNumberOfSets) {
@@ -73,6 +83,9 @@ void computeEvTime(const std::vector<eventTimeTrack>& tracks, const std::vector<
 
         int index = trkIndex[trackInSet[iset][itrk]];
         const eventTimeTrack& ctrack = tracks[trackInSet[iset][itrk]];
+        LOG(debug) << "Using hypothesis: " << hypo[itrk] << " tofSignal: " << ctrack.mSignal << " exp. time: " << ctrack.expTimes[hypo[itrk]] << " exp. sigma: " << ctrack.expSigma[hypo[itrk]];
+        LOG(debug) << "0= " << ctrack.expTimes[0] << " +- " << ctrack.expSigma[0] << " 1= " << ctrack.expTimes[1] << " +- " << ctrack.expSigma[1] << " 2= " << ctrack.expTimes[2] << " +- " << ctrack.expSigma[2];
+
         evtime.weights[index] = 1. / (ctrack.expSigma[hypo[itrk]] * ctrack.expSigma[hypo[itrk]]);
         evtime.tracktime[index] = ctrack.mSignal - ctrack.expTimes[hypo[itrk]];
       }
@@ -82,19 +95,21 @@ void computeEvTime(const std::vector<eventTimeTrack>& tracks, const std::vector<
   // do average among all tracks
   float finalTime = 0, allweights = 0;
   for (int i = 0; i < evtime.weights.size(); i++) {
-    if (evtime.weights[i] < 1E-6) {
+    if (evtime.weights[i] < weightLimit) {
       continue;
     }
     allweights += evtime.weights[i];
     finalTime += evtime.tracktime[i] * evtime.weights[i];
   }
 
-  if (allweights < 1E-6) {
+  if (allweights < weightLimit) {
+    LOG(debug) << "Skipping because allweights " << allweights << " are lower than " << weightLimit;
     return;
   }
 
   evtime.eventTime = finalTime / allweights;
   evtime.eventTimeError = sqrt(1. / allweights);
+  evtime.eventTimeMultiplicity = ntracks;
 }
 
 int getStartTimeInSet(const std::vector<eventTimeTrack>& tracks, std::vector<int>& trackInSet, unsigned long& bestComb)

Framework/Core/include/Framework/HistogramSpec.h

@@ -117,24 +117,7 @@ struct AxisSpec {
   }
 
   /// Function to make the axis logartitmic
-  void makeLogaritmic()
-  {
-    if (binEdges.size() > 2) {
-      LOG(FATAL) << "Cannot make a variabled bin width axis logaritmic";
-    }
-
-    const double min = binEdges[0];
-    const double max = binEdges[1];
-    const double logmin = std::log10(min);
-    const double logmax = std::log10(max);
-    const int nbins = nBins.value();
-    const double logdelta = (logmax - logmin) / (static_cast<double>(nbins));
-    const double log10 = std::log10(10.);
-    for (int i = 0; i < nbins + 1; i++) {
-      binEdges.push_back(std::exp(log10 * (logmin + i * logdelta)));
-    }
-    nBins = std::nullopt;
-  }
+  void makeLogaritmic();
 
   /// Data members
   std::optional<int> nBins{};         /// Number of bins (only used for fixed bin width axis)

Framework/Core/src/HistogramSpec.cxx

@@ -15,6 +15,29 @@
 namespace o2::framework
 {
 
+void AxisSpec::makeLogaritmic()
+{
+  if (binEdges.size() > 2) {
+    LOG(FATAL) << "Cannot make a variabled bin width axis logaritmic";
+  }
+
+  const double min = binEdges[0];
+  const double max = binEdges[1];
+  binEdges.clear();
+  const double logmin = std::log10(min);
+  const double logmax = std::log10(max);
+  const int nbins = nBins.value();
+  const double logdelta = (logmax - logmin) / (static_cast<double>(nbins));
+  const double log10 = std::log10(10.);
+  LOG(debug) << "Making a logaritmic binning from " << min << " to " << max << " with " << nbins << " bins";
+  for (int i = 0; i < nbins + 1; i++) {
+    const auto nextEdge = std::pow(10, logmin + i * logdelta);
+    LOG(debug) << i << "/" << nbins - 1 << ": " << nextEdge;
+    binEdges.push_back(nextEdge);
+  }
+  nBins = std::nullopt;
+}
+
 // main function for creating arbirtary histograms
 template <typename T>
 std::unique_ptr<T> HistFactory::createHist(const HistogramSpec& histSpec)