add run_reco_t0

AliceO2Group · sawenzel · Apr 1, 2019 · Jan 18, 2019 · Jan 22, 2019 · Jan 23, 2019
commit 15dc66a45710d68eafd1eaa522d751be4ad86352
@@ -29,7 +29,7 @@ CollisionTimeRecoTask::CollisionTimeRecoTask()
 void CollisionTimeRecoTask::Process(const Digit& digits,
                                     RecPoints& recPoints) const
 {
-  LOG(INFO) << "Running reconstruction on new event" << FairLogger::endl;
+  //  LOG(INFO) << "Running reconstruction on new event" << FairLogger::endl;
   recPoints.FillFromDigits(digits);
 }
 //________________________________________________________

@@ -21,14 +21,12 @@ ClassImp(o2::t0::RecPoints);
 
 void RecPoints::FillFromDigits(const o2::fit::Digit& digit)
 {
-  mTimeAmp.clear();
   mCollisionTime = {};
 
   Int_t ndigitsC = 0, ndigitsA = 0;
   constexpr Int_t nMCPsA = 4 * o2::t0::Geometry::NCellsA;
   constexpr Int_t nMCPsC = 4 * o2::t0::Geometry::NCellsC;
   constexpr Int_t nMCPs = nMCPsA + nMCPsC;
-  Float_t cfd[nMCPs] = {}, amp[nMCPs] = {};
   Float_t sideAtime = 0, sideCtime = 0;
 
   mBC = digit.getBC();
@@ -37,39 +35,28 @@ void RecPoints::FillFromDigits(const o2::fit::Digit& digit)
 
   Float_t BCEventTime = 12.5;
 
-  for (const auto& d : digit.getChDgData()) {
-    Int_t mcp = d.ChId;
-    cfd[mcp] = d.CFDTime - mEventTime - BCEventTime;
-    amp[mcp] = d.QTCAmpl;
-    mTimeAmp.push_back(ChannelData{ mcp, cfd[mcp], amp[mcp] });
-    //  LOG(DEBUG) << " mcp " << mcp<<" time "<< cfd[mcp]<<" amplitude "<< amp[mcp] << FairLogger::endl;
-  }
-
-  for (Int_t imcp = 0; imcp < nMCPsA; imcp++) {
-    if (cfd[imcp] > -2 && cfd[imcp] < 2) {
-      sideAtime += (cfd[imcp]);
-      ndigitsA++;
-    }
-  }
-  for (Int_t imcp = 0; imcp < nMCPsC; imcp++) {
-    if (cfd[imcp + nMCPsA] > 0) {
-      sideCtime += (cfd[imcp + nMCPsA]);
-      ndigitsC++;
+  mTimeAmp = digit.getChDgData();
+  for (auto& d : mTimeAmp) {
+    d.CFDTime -= mEventTime /*- BCEventTime*/;
+    if (abs(d.CFDTime - BCEventTime) < 2) {
+      if (d.ChId < nMCPsA) {
+        sideAtime += d.CFDTime;
+        ndigitsA++;
+      } else {
+        sideCtime += d.CFDTime;
+        ndigitsC++;
+      }
     }
   }
+
   if (ndigitsA > 0)
-    sideAtime = sideAtime / Float_t(ndigitsA);
+    mCollisionTime[1] = sideAtime / Float_t(ndigitsA);
+
   if (ndigitsC > 0)
-    sideCtime = sideCtime / Float_t(ndigitsC);
+    mCollisionTime[2] = sideCtime / Float_t(ndigitsC);
 
-  if (sideAtime > 0 && sideCtime > 0) {
-    mVertex = (sideAtime - sideCtime) / 2.;
-    mCollisionTime[0] = (sideAtime + sideCtime) / 2.;
-  }
-  if (sideAtime > 0) {
-    mCollisionTime[1] = sideAtime;
-  }
-  if (sideCtime > 0) {
-    mCollisionTime[2] = sideCtime;
+  if (ndigitsA > 0 && ndigitsC > 0) {
+    mVertex = (mCollisionTime[1] - mCollisionTime[2]) / 2.;
+    mCollisionTime[0] = (mCollisionTime[1] + mCollisionTime[2]) / 2.;
   }
 }
@@ -16,6 +16,7 @@
 
 #pragma link C++ class o2::t0::CollisionTimeRecoTask + ;
 #pragma link C++ class o2::t0::RecPoints + ;
+#pragma link C++ class std::vector < o2::t0::RecPoints > +;
 #pragma link C++ class o2::t0::Channel + ;
 
 #endif
@@ -0,0 +1,72 @@
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include <Rtypes.h>
+#include <TString.h>
+#include <TStopwatch.h>
+#include <TGeoManager.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <iostream>
+
+#include "FairLogger.h"
+
+#include "T0Reconstruction/CollisionTimeRecoTask.h"
+#include "T0Reconstruction/RecPoints.h"
+#include "FITBase/Digit.h"
+#endif
+
+void run_reco_t0(std::string inpudDig = "t0digits.root",
+                 std::string outName = "o2reco_t0.root",
+                 std::string inputGRP = "o2sim_grp.root")
+{
+
+  // Initialize logger
+  FairLogger* logger = FairLogger::GetLogger();
+  logger->SetLogVerbosityLevel("LOW");
+  logger->SetLogScreenLevel("DEBUG");
+
+  // Setup timer
+  TStopwatch timer;
+
+  TFile* fdig = TFile::Open("t0digits.root");
+  std::cout << " Open digits file " << std::endl;
+  TTree* digTree = (TTree*)fdig->Get("o2sim");
+  digTree->Print();
+  std::vector<o2::fit::Digit>* digits = new std::vector<o2::fit::Digit>;
+  digTree->SetBranchAddress("T0Digit", &digits);
+  Int_t nevD = digTree->GetEntries(); // digits in cont. readout may be grouped as few events per entry
+  std::cout << "Found " << nevD << " events with digits " << std::endl;
+
+  //std::vector<o2::fit::RecPoints> recPoints;
+  std::vector<o2::t0::RecPoints>* recPointsP = new std::vector<o2::t0::RecPoints>;
+
+  TFile outFile(outName.c_str(), "recreate");
+  TTree outTree("o2sim", "T0RecPoints");
+  outTree.Branch("T0Cluster", &recPointsP);
+
+  o2::t0::CollisionTimeRecoTask recoFIT;
+  timer.Start();
+  for (int iEv = 0; iEv < nevD; iEv++) {
+    digTree->GetEntry(iEv);
+    recPointsP->resize(digits->size());
+    for (size_t collID = 0; collID < digits->size(); ++collID)
+      recoFIT.Process((*digits)[collID], (*recPointsP)[collID]);
+    outTree.Fill();
+  }
+
+  outTree.Print();
+  outFile.ls();
+  timer.Stop();
+  outTree.Write();
+  Double_t rtime = timer.RealTime();
+  Double_t ctime = timer.CpuTime();
+
+  Float_t cpuUsage = ctime / rtime;
+  std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+  std::cout << cpuUsage;
+  std::cout << "</DartMeasurement>" << std::endl;
+  std::cout << "Macro finished succesfully." << std::endl;
+
+  std::cout << "Real time " << rtime << " s, CPU time " << ctime
+            << "s" << std::endl
+            << std::endl;
+}