add first tof start time prototype

Author: njacazio

Detectors/TOF/prototyping/CMakeLists.txt

@@ -15,7 +15,7 @@ o2_add_test_root_macro(checkRotation.C
 
 o2_add_test_root_macro(convertTreeTo02object.C
                        PUBLIC_LINK_LIBRARIES O2::ReconstructionDataFormats
-		       			     O2::DataFormatsTOF 
+                                             O2::DataFormatsTOF
                        LABELS tof)
 
 o2_add_test_root_macro(drawTOFgeometry.C
@@ -43,14 +43,19 @@ o2_add_test_root_macro(findTOFclusterFromLabel.C
 
 o2_add_test_root_macro(checkTS.C
                        PUBLIC_LINK_LIBRARIES O2::DataFormatsTOF
-		       LABELS tof)
+                       LABELS tof)
 
 o2_add_test_root_macro(testCCDB.C
                        PUBLIC_LINK_LIBRARIES O2::DataFormatsTOF
-		       			     O2::TOFCalibration
-		       LABELS tof)
+                                             O2::TOFCalibration
+                       LABELS tof)
 
 o2_add_test_root_macro(checkChannelMapping.C
                        PUBLIC_LINK_LIBRARIES O2::TOFBase
                        LABELS tof)
 
+
+o2_add_test_root_macro(macroEvTime.C
+                       PUBLIC_LINK_LIBRARIES O2::TOFBase
+                                             O2::TOFReconstruction
+                       LABELS tof)

Detectors/TOF/prototyping/macroEvTime.C

@@ -0,0 +1,61 @@
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include "TOFReconstruction/EventTimeMaker.h"
+#endif
+
+using namespace o2::tof;
+void macroEvTime(bool removebias = true)
+{
+
+  TH1F* htest = new TH1F("htest", ";N#sigma (t0_{reco} - t0_{true})", 100, -10, 10);
+  TH1F* htest2 = new TH1F("htest2", ";t0_{reco} - t0_{true} ps", 100, -100, 100);
+
+  TH1F* hkaon = new TH1F("hkaon", ";N#sigma^{K}", 100, -10, 10);
+  TH1F* hkaonT = new TH1F("hkaonT", ";N#sigma^{K}", 100, -10, 10);
+  hkaonT->SetLineColor(2);
+
+  std::vector<eventTimeTrackTest> tracks;
+  for (int i = 0; i < 1000; i++) {
+    tracks.clear();
+    generateEvTimeTracks(tracks, 5);
+    auto evtime = evTimeMaker<std::vector<eventTimeTrackTest>, eventTimeTrackTest, filterDummy>(tracks);
+    //    Printf("Ev time %f +-%f", evtime.eventTime, evtime.eventTimeError);
+    htest->Fill(evtime.eventTime / evtime.eventTimeError);
+    htest2->Fill(evtime.eventTime);
+
+    int nt = 0;
+    for (auto track : tracks) {
+
+      if (track.p() > 2) {
+        nt++;
+        continue;
+      }
+
+      float et = evtime.eventTime;
+      float erret = evtime.eventTimeError;
+
+      if (removebias) {
+        float sumw = 1. / erret / erret;
+        et *= sumw;
+        et -= evtime.weights[nt] * evtime.tracktime[nt];
+        sumw -= evtime.weights[nt];
+        et /= sumw;
+        erret = sqrt(1. / sumw);
+      }
+
+      nt++;
+
+      float res = sqrt(100 * 100 + erret * erret);
+      hkaon->Fill((track.tofSignal() - track.tofExpTimeKa() - et) / res);
+
+      if (track.masshypo() == 1) {
+        hkaonT->Fill((track.tofSignal() - track.tofExpTimeKa() - et) / res);
+      }
+    }
+  }
+  htest->Draw();
+  new TCanvas;
+  htest2->Draw();
+  new TCanvas;
+  hkaon->Draw();
+  hkaonT->Draw("SAME");
+}

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

@@ -15,6 +15,9 @@
 #ifndef ALICEO2_TOF_EVENTTIMEMAKER_H
 #define ALICEO2_TOF_EVENTTIMEMAKER_H
 
+#include "TRandom.h"
+#include "TMath.h"
+
 namespace o2
 {
 
@@ -25,18 +28,93 @@ 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
+};
+
+struct eventTimeTrack {
+  eventTimeTrack() = default;
+  eventTimeTrack(float tof, float expt[3], float expsigma[3]) : mSignal(tof)
+  {
+    for (int i = 0; i < 3; i++) {
+      expTimes[i] = expt[i];
+      expSigma[i] = expsigma[i];
+    }
+  }
+  float tofSignal() const { return mSignal; }
+  float tofExpTimePi() const { return expTimes[0]; }
+  float tofExpTimeKa() const { return expTimes[1]; }
+  float tofExpTimePr() const { return expTimes[2]; }
+  float tofExpSigmaPi() const { return expSigma[0]; }
+  float tofExpSigmaKa() const { return expSigma[1]; }
+  float tofExpSigmaPr() const { return expSigma[2]; }
+  float mSignal = 0.f;
+  float expTimes[3] = {0.f, 0.f, 0.f};
+  float expSigma[3] = {999.f, 999.f, 999.f};
+};
+
+struct eventTimeTrackTest : eventTimeTrack {
+  float tofChi2() const { return mTOFChi2; }
+  float pt() const { return mPt; }
+  float p() const { return mP; }
+  float length() const { return mLength; }
+  int masshypo() const { return mHypo; }
+  float mTOFChi2 = -1.f;
+  float mPt = 0.f;
+  float mP = 0.f;
+  float mLength = 0.f;
+  int mHypo = 0;
 };
 
-template <typename trackContainer>
-eventTimeContainer evTimeMaker(const trackContainer& tracks)
+void generateEvTimeTracks(std::vector<eventTimeTrackTest>& tracks, int ntracks, float evTime = 0.f);
+
+template <typename trackType>
+bool filterDummy(const trackType& tr)
 {
+  return (tr.tofChi2() >= 0 && tr.mP < 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 = 6. /* 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();
+
+  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.;
+
+  // Qui facciamo un pool di tracce buone per calcolare il T0
   for (auto track : tracks) {
-    track.tofSignal();
+    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();
+      trkWork.emplace_back(track.tofSignal(), expt, expsigma);
+      trkIndex.push_back(result.weights.size());
+    }
+    result.weights.push_back(0.);
+    result.tracktime.push_back(0.);
   }
-  return eventTimeContainer{0.f, 0.f};
+  computeEvTime(trkWork, trkIndex, result);
+  return result;
 }
 
 } // namespace tof
 } // namespace o2
 
-#endif /* ALICEO2_TOF_EVENTTIMEMAKER_H */
+#endif /* ALICEO2_TOF_EVENTTIMEMAKER_H */