diff --git a/reco/detectors/trd/CbmTrdClusterFinder.h b/reco/detectors/trd/CbmTrdClusterFinder.h index 7ae8d2a71e8a4a12aa5ce6bb23df3438d4ebfbfd..99a48d1ee8818034d4b5a55b7c3a1fcb997f9f14 100644 --- a/reco/detectors/trd/CbmTrdClusterFinder.h +++ b/reco/detectors/trd/CbmTrdClusterFinder.h @@ -58,17 +58,18 @@ class CbmTrdClusterFinder : public FairTask { friend class CbmTrdModuleRec2D; public: - enum CbmTrdRecDef + enum ECbmTrdRecDef { - kTime = 0, ///< select Time based/Event by event reconstruction - kMultiHit, ///< multi hit detection - kRowMerger, ///< merge clusters over neighbour rows - kNeighbourCol, ///< use neighbour trigger; column wise - kNeighbourRow, ///< use neighbour trigger; row wise - kDumpClusters, ///< write clustered digis to output - kFASP, ///< use FASP ASIC for triangular pad plane geometry - kOnlyEventDigis, ///< use only digis connected to a CbmEvent - kDebugStatements ///< print out debug statements with LOG(info) + kTime = 0, ///< select Time based/Event by event reconstruction + kMultiHit, ///< multi hit detection + kRowMerger, ///< merge clusters over neighbor rows + kNeighbourCol, ///< use neighbor trigger; column wise + kNeighbourRow, ///< use neighbor trigger; row wise + kDumpClusters, ///< write clustered digis to output + kFASP, ///< use FASP ASIC for triangular pad plane geometry + kOnlyEventDigis, ///< use only digis connected to a CbmEvent + kDebugStatements, ///< print out debug statements with LOG(info) + kUseRecoHelpers ///< use graph helpers in hit reco for improving performance }; /** @@ -89,11 +90,12 @@ public: static Bool_t HasRowMerger() { return TESTBIT(fgConfig, kRowMerger); } static Bool_t IsTimeBased() { return TESTBIT(fgConfig, kTime); } static Bool_t DoDebugPrintouts() { return TESTBIT(fgConfig, kDebugStatements); } + static Bool_t UseRecoHelpers() { return TESTBIT(fgConfig, kUseRecoHelpers); } /** * @brief If true only digis connected ro a CbmEvent are processed * Per default this is activated on construction, such that the user can - * turn it off before Init(), where it will also be autmatically + * turn it off before Init(), where it will also be automatically * deactivated if no CbmEvent branch is found. * @return Bool_t */ @@ -129,12 +131,19 @@ public: { set ? SETBIT(fgConfig, kDebugStatements) : CLRBIT(fgConfig, kDebugStatements); } + /** + * @brief Steer usage of TGraph support in TRD-2D for improved performance. The cost in terms of CPU has to be weighted against performance increase + */ + static void SetUseRecoHelpers(Bool_t set = kTRUE) + { + set ? SETBIT(fgConfig, kUseRecoHelpers) : CLRBIT(fgConfig, kUseRecoHelpers); + } /** * @brief Set the Use Only Event Digis * Per default this is activated on construction, such that the user can - * turn it off before Init(), where it will also be autmatically + * turn it off before Init(), where it will also be automatically * deactivated if no CbmEvent branch is found. * @param set */ diff --git a/reco/detectors/trd/CbmTrdHitProducer.cxx b/reco/detectors/trd/CbmTrdHitProducer.cxx index 79459969e93bdacf8fead09d9f5a3d437bfa3b5e..f7a2e68caeb0af2987c4ea0fa8806773074a4383 100644 --- a/reco/detectors/trd/CbmTrdHitProducer.cxx +++ b/reco/detectors/trd/CbmTrdHitProducer.cxx @@ -98,7 +98,10 @@ CbmTrdModuleRec* CbmTrdHitProducer::AddModule(Int_t address, TGeoPhysicalNode* n << moduleAddress << "] ly[" << lyId << "] det[" << moduleAddress << "]"; CbmTrdModuleRec* module(nullptr); - if (moduleType >= 9) { module = fModules[address] = new CbmTrdModuleRec2D(address); } + if (moduleType >= 9) { + module = fModules[address] = new CbmTrdModuleRec2D(address); + ((CbmTrdModuleRec2D*) module)->SetUseHelpers(CbmTrdClusterFinder::UseRecoHelpers()); + } else { module = fModules[address] = new CbmTrdModuleRecR(address); } @@ -339,7 +342,8 @@ void CbmTrdHitProducer::Exec(Option_t*) logOut << fixed << setw(8) << setprecision(1) << right << timerTs.RealTime() * 1000. << " ms] "; logOut << "TS " << fNrTs; if (CbmTrdClusterFinder::UseOnlyEventDigis()) logOut << ", events " << nEvents; - logOut << ", clusters " << nClusters << ", hits " << fNrHitsCall; + logOut << ", clusters " << nClusters << ", hits " << fNrHitsCall << ", time/1k-hit " << setprecision(4) + << timerTs.RealTime() * 1e6 / fNrHitsCall << " [ms]"; LOG(info) << logOut.str(); fNrTs++; } diff --git a/reco/detectors/trd/CbmTrdModuleRec2D.cxx b/reco/detectors/trd/CbmTrdModuleRec2D.cxx index de121eed68fa518d1f1c0269bad6cf24bfd5ecf1..5ae7c93b1bdbac1ff20bf9fd02c783e05108e07b 100644 --- a/reco/detectors/trd/CbmTrdModuleRec2D.cxx +++ b/reco/detectors/trd/CbmTrdModuleRec2D.cxx @@ -642,74 +642,88 @@ Bool_t CbmTrdModuleRec2D::BuildHit(CbmTrdHit* h) LocalToMaster(local, global); // COMPUTE TIME + Double_t t_avg(0.), e_avg(0.); for (Int_t idx(1); idx <= n0; idx++) { Double_t dtFEE = fgDT[0] * (vs[idx] - fgDT[1]) * (vs[idx] - fgDT[1]) / CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kFASP); if (vxe[idx] > 0) vx[idx] += dy / fDigiPar->GetPadSizeY(0); - fgT->SetPoint(idx - 1, vx[idx], vt[idx] - dtFEE); + if (fgT != nullptr) fgT->SetPoint(idx - 1, vx[idx], vt[idx] - dtFEE); + else + t_avg += (vt[idx] - dtFEE); } - Double_t xc = vx[n0 + 2]; - for (Int_t ip(n0); ip < fgT->GetN(); ip++) { - fgT->SetPoint(ip, xc, 0); - xc += 0.5; + Double_t xc(0.); + if (fgT != nullptr) { + xc = vx[n0 + 2]; + for (Int_t ip(n0); ip < fgT->GetN(); ip++) { + fgT->SetPoint(ip, xc, 0); + xc += 0.5; + } } Double_t time(-21.), tdrift(100); // should depend on Ua - if (n0 > 1 && (fgT->Fit("pol1", "QC", "goff") == 0)) { - TF1* f = fgT->GetFunction("pol1"); - time = f->GetParameter(0) - fgDT[2]; - if (TMath::IsNaN(time)) time = -21; - //dtime += TMath::Abs(f->GetParameter(1)*(vx[n0+1] - vx[1])); + if (n0 > 1) { + if (fgT != nullptr && (fgT->Fit("pol1", "QC", "goff") == 0)) { + TF1* f = fgT->GetFunction("pol1"); + time = f->GetParameter(0) - fgDT[2]; + if (TMath::IsNaN(time)) time = -21; + //dtime += TMath::Abs(f->GetParameter(1)*(vx[n0+1] - vx[1])); + } + else + time = t_avg / n0; } - // COMPUTE ENERGY for (UInt_t idx(0); idx < vs.size(); idx++) { - if (vxe[idx] > 0) { - fgEdep->SetPoint(idx, vx[idx] + dy / fDigiPar->GetPadSizeY(0), vs[idx]); - fgEdep->SetPointError(idx, vxe[idx], vse[idx]); - } - else { - fgEdep->SetPoint(idx, vx[idx], vs[idx]); + if (fgEdep != nullptr) { + double x_offset = dy / fDigiPar->GetPadSizeY(0), xp = vx[idx] + (vxe[idx] > 0 ? x_offset : 0); + fgEdep->SetPoint(idx, xp, vs[idx]); fgEdep->SetPointError(idx, vxe[idx], vse[idx]); } + else + e_avg += vs[idx]; } - xc = vx[n0 + 2]; - for (Int_t ip(vs.size()); ip < fgEdep->GetN(); ip++) { - //fgEdep->RemovePoint(ip); - xc += 0.5; - fgEdep->SetPoint(ip, xc, 0); - fgEdep->SetPointError(ip, 0., 300); + if (fgEdep != nullptr) { + xc = vx[n0 + 2]; + for (Int_t ip(vs.size()); ip < fgEdep->GetN(); ip++) { + //fgEdep->RemovePoint(ip); + xc += 0.5; + fgEdep->SetPoint(ip, xc, 0); + fgEdep->SetPointError(ip, 0., 300); + } + //if(CWRITE()) fgEdep->Print(); } - //if(CWRITE()) fgEdep->Print(); Double_t e(0.), xlo(*vx.begin()), xhi(*vx.rbegin()); - fgPRF->SetParameter(0, vs[viM]); - fgPRF->FixParameter(1, dx / fDigiPar->GetPadSizeX(0)); - fgPRF->SetParameter(2, 0.65); - fgPRF->SetParLimits(2, 0.45, 10.5); - fgEdep->Fit(fgPRF, "QBN", "goff", xlo - 0.5, xhi + 0.5); - if (!fgPRF->GetNDF()) return false; - //chi = fgPRF->GetChisquare()/fgPRF->GetNDF(); - e = fgPRF->Integral(xlo - 0.5, xhi + 0.5); - - // apply MC correction - Float_t gain0 = 210.21387; //(XeCO2 @ 1900V) - // Float_t grel[3] = {1., 0.98547803, 0.93164071}, - // goff[3][3] = { - // {0.05714, -0.09, -0.09}, - // {0., -0.14742, -0.14742}, - // {0., -0.29, -0.393} - // }; - // Int_t ian=0; - // if(TMath::Abs(dy)<=0.3) ian=0; - // else if(TMath::Abs(dy)<=0.6) ian=1; - // else if(TMath::Abs(dy)<=0.9) ian=2; - // Int_t isize=0; - // if(n0<=3) isize=0; - // else if(n0<=4) isize=1; - // else isize=2; - Float_t gain = gain0; //*grel[ian]; - e /= gain; // apply effective gain - //e+=goff[ian][isize]; // apply missing energy offset + if (fgEdep != nullptr) { + fgPRF->SetParameter(0, vs[viM]); + fgPRF->FixParameter(1, dx / fDigiPar->GetPadSizeX(0)); + fgPRF->SetParameter(2, 0.65); + fgPRF->SetParLimits(2, 0.45, 10.5); + fgEdep->Fit(fgPRF, "QBN", "goff", xlo - 0.5, xhi + 0.5); + if (!fgPRF->GetNDF()) return false; + //chi = fgPRF->GetChisquare()/fgPRF->GetNDF(); + e = fgPRF->Integral(xlo - 0.5, xhi + 0.5); + + // apply MC correction + Float_t gain0 = 210.21387; //(XeCO2 @ 1900V) + // Float_t grel[3] = {1., 0.98547803, 0.93164071}, + // goff[3][3] = { + // {0.05714, -0.09, -0.09}, + // {0., -0.14742, -0.14742}, + // {0., -0.29, -0.393} + // }; + // Int_t ian=0; + // if(TMath::Abs(dy)<=0.3) ian=0; + // else if(TMath::Abs(dy)<=0.6) ian=1; + // else if(TMath::Abs(dy)<=0.9) ian=2; + // Int_t isize=0; + // if(n0<=3) isize=0; + // else if(n0<=4) isize=1; + // else isize=2; + Float_t gain = gain0; //*grel[ian]; + e /= gain; // apply effective gain + //e+=goff[ian][isize]; // apply missing energy offset + } + else + e = e_avg; h->SetX(global[0]); h->SetY(global[1]); @@ -746,7 +760,7 @@ CbmTrdHit* CbmTrdModuleRec2D::MakeHit(Int_t ic, const CbmTrdCluster* cl, std::ve } //printf("%s (%s)\n", GetName(), GetTitle()); Config(1,0); - if (!fgEdep) { // first use initialization of PRF helppers + if (CHELPERS() && fgEdep == nullptr) { // first use initialization of PRF helppers LOG(info) << GetName() << "::MakeHit: Init static helpers. "; fgEdep = new TGraphErrors; fgEdep->SetLineColor(kBlue); diff --git a/reco/detectors/trd/CbmTrdModuleRec2D.h b/reco/detectors/trd/CbmTrdModuleRec2D.h index c058c545329669bdb14aa6e1f9c6d58e71fe7dba..f74f0e9cb45079d3262c930e947726af0d35640e 100644 --- a/reco/detectors/trd/CbmTrdModuleRec2D.h +++ b/reco/detectors/trd/CbmTrdModuleRec2D.h @@ -35,11 +35,11 @@ class TF1; **/ class CbmTrdModuleRec2D : public CbmTrdModuleRec { public: - enum CbmTrdModuleRec2Dconfig + enum ECbmTrdModuleRec2D { - kVerbose = 0 // steer verbosity on/off - , - kDraw // steer graphic representation on/off + kVerbose = 0, ///< steer verbosity on/off + kDraw = 1, ///< steer graphic representation on/off + kHelpers = 2 ///< use helper graph for time and energy estimation }; /** \brief Default constructor.*/ @@ -118,13 +118,23 @@ public: */ Int_t GetHitRcClass(Int_t a0) const; + /** @brief Steer usage of helper graphs for computing time and energy per hit. + * A cost wrt the additional performance vs. CPU has to be performed. + */ + void SetUseHelpers(bool use = true) + { + use ? SETBIT(fConfigMap, ECbmTrdModuleRec2D::kHelpers) : CLRBIT(fConfigMap, ECbmTrdModuleRec2D::kHelpers); + } + protected: private: CbmTrdModuleRec2D(const CbmTrdModuleRec2D& ref); const CbmTrdModuleRec2D& operator=(const CbmTrdModuleRec2D& ref); - Bool_t CDRAW() const { return TESTBIT(fConfigMap, kDraw); } - Bool_t CWRITE() const { return TESTBIT(fConfigMap, kVerbose); } + Bool_t CDRAW() const { return TESTBIT(fConfigMap, ECbmTrdModuleRec2D::kDraw); } + Bool_t CWRITE() const { return TESTBIT(fConfigMap, ECbmTrdModuleRec2D::kVerbose); } + Bool_t CHELPERS() const { return TESTBIT(fConfigMap, ECbmTrdModuleRec2D::kHelpers); } + /** \brief Implement cuts for hit convolution definition * \param[in] h hit to be analysed. * \return TRUE if double cluster @@ -217,13 +227,13 @@ private: void CbmTrdModuleRec2D::Config(Bool_t vb, Bool_t dw) { - if (vb) SETBIT(fConfigMap, kVerbose); + if (vb) SETBIT(fConfigMap, ECbmTrdModuleRec2D::kVerbose); else - CLRBIT(fConfigMap, kVerbose); + CLRBIT(fConfigMap, ECbmTrdModuleRec2D::kVerbose); printf("CbmTrdModuleRec2D::Verbose[%c]\n", CWRITE() ? 'y' : 'n'); - if (dw) SETBIT(fConfigMap, kDraw); + if (dw) SETBIT(fConfigMap, ECbmTrdModuleRec2D::kDraw); else - CLRBIT(fConfigMap, kDraw); + CLRBIT(fConfigMap, ECbmTrdModuleRec2D::kDraw); printf("CbmTrdModuleRec2D::Draw[%c]\n", CDRAW() ? 'y' : 'n'); }