diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx index bdf45e8fef6e22f7ffc5e00ce45a1bf34108c49c..ffa17d451ab99ff3375748e64136a2da5d559ed7 100644 --- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx +++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx @@ -61,9 +61,7 @@ Bool_t CbmMcbmCheckTimingAlgo::Init() /// Try to get the 2021 Event header which is containing the Timeslice info /// If not present, we have "old" data and will simply catch it with the nullptr value - auto eh = FairRun::Instance()->GetEventHeader(); - LOG(info) << "CbmMcbmCheckTimingAlgo => EventHeader ptr: " << eh << " (" << eh->IsA()->GetName() << ")"; - fCbmTsEventHeader = dynamic_cast<CbmTsEventHeader*>(eh); + fCbmTsEventHeader = FairRootManager::Instance()->InitObjectAs<CbmTsEventHeader const*>("EventHeader."); if (nullptr != fCbmTsEventHeader) { LOG(info) << "CbmMcbmCheckTimingAlgo => Using the index from the TS Event header for Evo plots"; } @@ -105,7 +103,7 @@ void CbmMcbmCheckTimingAlgo::CreateHistos() for( std::vector< CheckTimingDetector >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) { - for (uint i(0); i < (*det).uNdiv; i++) { + for (uint i(0); i < (*det).uNviews; i++) { fvhDetSelfDiff[(*det).detId].push_back( new TH1D(Form("h%sSelfDiff%d", (*det).sName.data(), i), Form("time difference between consecutivs %s (%s) Digis;time diff [ns];Counts", (*det).sName.data(), @@ -182,8 +180,6 @@ void CbmMcbmCheckTimingAlgo::ProcessTs() LOG(info) << "executing TS " << fuNbTs << " index [" << (fCbmTsEventHeader != nullptr ? fCbmTsEventHeader->GetTsIndex() : -1) << "]"; - fuNbTrdCondWarn = 0; - switch (fRefDet.detId) { case ECbmModuleId::kSts: { CheckInterSystemOffset<CbmStsDigi>(); @@ -358,61 +354,54 @@ void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const UInt_t uFirstDigiInWin = fvDets[uDetIdx].iPrevRefFirstDigi; - std::vector<Double_t> vSelDiff; + bool exit = false; + std::vector<Double_t> vSelDiff = {}; + std::vector<std::tuple<double, double, uint>> vDigiInfo = {}; // digi info containing (time, charge, address) for (UInt_t uDigiIdx = fvDets[uDetIdx].iPrevRefFirstDigi; uDigiIdx < uNbDigis; ++uDigiIdx) { - Double_t dTime = 0; - Double_t dCharge = 0; - UInt_t uAddress = 0; - const void* digi(nullptr); - if (ECbmModuleId::kT0 == edetId) { - dTime = fpT0DigiVec->at(uDigiIdx).GetTime(); - dCharge = fpT0DigiVec->at(uDigiIdx).GetCharge(); - } // if( ECbmModuleId::kT0 == fRefDet.detId ) - else { - digi = fDigiMan->Get<Digi>(uDigiIdx); - dTime = fDigiMan->Get<Digi>(uDigiIdx)->GetTime(); - dCharge = fDigiMan->Get<Digi>(uDigiIdx)->GetCharge(); - uAddress = fDigiMan->Get<Digi>(uDigiIdx)->GetAddress(); - } // else of if( ECbmModuleId::kT0 == fRefDet.detId ) - - /// Fill self correlation histo while avoiding double counting due to - /// the "smart looping" - if (fvDets[uDetIdx].dPrevTime <= dTime) { - vSelDiff.push_back(dTime - fvDets[uDetIdx].dPrevTime); - fvDets[uDetIdx].dPrevTime = dTime; - } // if( fvDets[ uDetIdx ].dPrevTime < dTime ) - - Double_t dDiffTime = dTime - dRefTime; - - if (dDiffTime < fvDets[uDetIdx].dTimeRangeBeg) { - ++uFirstDigiInWin; // Update Index of first digi in Win to next digi - continue; // not yet in interesting range - } // if (diffTime > offsetRange) - if (fvDets[uDetIdx].dTimeRangeEnd < dDiffTime) { - /// already past interesting range - break; - } // if( fvDets[ uDetIdx ].dTimeRangeEnd < dDiffTime ) - - /// Charge cut if defined! - if (fvDets[uDetIdx].uChargeCutMin != fvDets[uDetIdx].uChargeCutMax) { - if (fvDets[uDetIdx].uChargeCutMin < fvDets[uDetIdx].uChargeCutMax) { - /// Cut charges between Min and Max to reject pulser - if (fvDets[uDetIdx].uChargeCutMin < dCharge && dCharge < fvDets[uDetIdx].uChargeCutMax) { - continue; - } // if( fvDets[ uDetIdx ].uChargeCutMin < dCharge && dCharge < fvDets[ uDetIdx ].uChargeCutMax ) - } // if( fvDets[ uDetIdx ].uChargeCutMin < fvDets[ uDetIdx ].uChargeCutMax ) - else { - /// Select charges between Max and Min to select pulser (Min and Max swapped!!) - if (fvDets[uDetIdx].uChargeCutMin < dCharge || dCharge < fvDets[uDetIdx].uChargeCutMax) { - continue; - } // if( fvDets[ uDetIdx ].uChargeCutMin < dCharge || dCharge < fvDets[ uDetIdx ].uChargeCutMax ) - /// Psd Pulser selection - if (ECbmModuleId::kPsd == fvDets[uDetIdx].detId && CbmPsdAddress::GetSectionId(uAddress) != 10) continue; - } // else of if( fvDets[ uDetIdx ].uChargeCutMin < fvDets[ uDetIdx ].uChargeCutMax ) - } // if( fvDets[ uDetIdx ].uChargeCutMin != fvDets[ uDetIdx ].uChargeCutMax ) - - for (uint uid(0); uid < fvDets[uDetIdx].uNdiv; uid++) { - if (digi && !CheckCondition(uDetIdx, uid, digi)) continue; + if (!GetDigiInfo<Digi>(uDigiIdx, &vDigiInfo, edetId)) continue; + for (auto dInfo : vDigiInfo) { + Double_t dTime = std::get<0>(dInfo); + Double_t dCharge = std::get<1>(dInfo); + UInt_t uAddress = std::get<2>(dInfo); + + /// Fill self correlation histo while avoiding double counting due to + /// the "smart looping" + if (fvDets[uDetIdx].dPrevTime <= dTime) { + vSelDiff.push_back(dTime - fvDets[uDetIdx].dPrevTime); + fvDets[uDetIdx].dPrevTime = dTime; + } // if( fvDets[ uDetIdx ].dPrevTime < dTime ) + Double_t dDiffTime = dTime - dRefTime; + if (dDiffTime < fvDets[uDetIdx].dTimeRangeBeg) { + ++uFirstDigiInWin; // Update Index of first digi in Win to next digi + //continue; // not yet in interesting range + break; // not yet in interesting range + } // if (diffTime > offsetRange) + if (fvDets[uDetIdx].dTimeRangeEnd < dDiffTime) { + exit = true; + /// already past interesting range + break; + } // if( fvDets[ uDetIdx ].dTimeRangeEnd < dDiffTime ) + + /// Charge cut if defined! + if (fvDets[uDetIdx].uChargeCutMin != fvDets[uDetIdx].uChargeCutMax) { + if (fvDets[uDetIdx].uChargeCutMin < fvDets[uDetIdx].uChargeCutMax) { + /// Cut charges between Min and Max to reject pulser + if (fvDets[uDetIdx].uChargeCutMin < dCharge && dCharge < fvDets[uDetIdx].uChargeCutMax) { + continue; + } // if( fvDets[ uDetIdx ].uChargeCutMin < dCharge && dCharge < fvDets[ uDetIdx ].uChargeCutMax ) + } // if( fvDets[ uDetIdx ].uChargeCutMin < fvDets[ uDetIdx ].uChargeCutMax ) + else { + /// Select charges between Max and Min to select pulser (Min and Max swapped!!) + if (fvDets[uDetIdx].uChargeCutMin < dCharge || dCharge < fvDets[uDetIdx].uChargeCutMax) { + continue; + } // if( fvDets[ uDetIdx ].uChargeCutMin < dCharge || dCharge < fvDets[ uDetIdx ].uChargeCutMax ) + /// Psd Pulser selection + if (ECbmModuleId::kPsd == fvDets[uDetIdx].detId && CbmPsdAddress::GetSectionId(uAddress) != 10) continue; + } // else of if( fvDets[ uDetIdx ].uChargeCutMin < fvDets[ uDetIdx ].uChargeCutMax ) + } // if( fvDets[ uDetIdx ].uChargeCutMin != fvDets[ uDetIdx ].uChargeCutMax ) + + int uid = GetViewId<Digi>(fvDets[uDetIdx], dInfo); + if (uid < 0 || uint(uid) >= fvDets[uDetIdx].uNviews) continue; /// Fill histos for (auto dt : vSelDiff) @@ -430,65 +419,88 @@ void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const fvhDetToRefDiffEvo[edetId][uid]->Fill(fCbmTsEventHeader->GetTsIndex(), dDiffTime); fvhDetToRefDiffEvoLong[edetId][uid]->Fill(fCbmTsEventHeader->GetTsIndex(), dDiffTime); } - break; - } + } // for (auto dInfo : vDigiInfo) { + if (exit) break; } // for( UInt_t uDigiIdx = fvDets[ uDetIdx ].iPrevRefFirstDigi; uDigiIdx < uNbDigis; ++uDigiIdx ) /// Store earliest possible starting index for next reference digi (time sorted!) fvDets[uDetIdx].iPrevRefFirstDigi = uFirstDigiInWin; } -// ---- CheckCondition ------------------------------------------------ -bool CbmMcbmCheckTimingAlgo::CheckCondition(UInt_t uDetId, UInt_t uCond, const void* digi) -{ - CheckTimingDetector uDet = fvDets[uDetId]; - if (uCond >= uDet.uNdiv) return false; - switch (uDet.detId) { - case ECbmModuleId::kSts: return CheckCondition(&uDet, uCond, static_cast<const CbmStsDigi*>(digi)); - case ECbmModuleId::kMuch: return CheckCondition(&uDet, uCond, static_cast<const CbmMuchDigi*>(digi)); - case ECbmModuleId::kTrd: return CheckCondition(&uDet, uCond, static_cast<const CbmTrdDigi*>(digi)); - case ECbmModuleId::kTof: return CheckCondition(&uDet, uCond, static_cast<const CbmTofDigi*>(digi)); - case ECbmModuleId::kRich: return CheckCondition(&uDet, uCond, static_cast<const CbmRichDigi*>(digi)); - case ECbmModuleId::kPsd: return CheckCondition(&uDet, uCond, static_cast<const CbmPsdDigi*>(digi)); - default: break; - } - return true; -} -bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmStsDigi* /*digi*/) +// ---- GetDigiTime ------------------------------------------------------ +template<class Digi> +uint CbmMcbmCheckTimingAlgo::GetDigiInfo(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec, ECbmModuleId) { - return true; + const Digi* digi = fDigiMan->Get<Digi>(uDigi); + vec->clear(); + if (digi == nullptr) return 0; + vec->push_back(std::make_tuple(digi->GetTime(), digi->GetCharge(), digi->GetAddress())); + return 1; } -bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmMuchDigi* /*digi*/) +template<> +uint CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmTofDigi>(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec, + ECbmModuleId detId) { - return true; + /** Template specialization for ToF in order to distinguish T0 digis. + */ + vec->clear(); + const CbmTofDigi* pDigi(nullptr); + if (detId == ECbmModuleId::kT0) pDigi = &fpT0DigiVec->at(uDigi); + else + pDigi = fDigiMan->Get<CbmTofDigi>(uDigi); + + if (pDigi == nullptr) return 0; + vec->push_back(std::make_tuple(pDigi->GetTime(), pDigi->GetCharge(), pDigi->GetAddress())); + return 1; } -bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmTrdDigi* digi) +template<> +uint CbmMcbmCheckTimingAlgo::GetDigiInfo<CbmTrdDigi>(UInt_t uDigi, std::vector<std::tuple<double, double, uint>>* vec, + ECbmModuleId) { - try { - int moduleId = std::stoi(det->vName[uCond]); - if (digi->GetAddressModule() != moduleId) return false; + /** Template specialization for TRD in order to distinguish SPADIC and FASP digis. + */ + vec->clear(); + const CbmTrdDigi* trdDigi = fDigiMan->Get<CbmTrdDigi>(uDigi); + if (trdDigi == nullptr) return 0; + + const double dt = trdDigi->GetTime(); + + if (trdDigi->GetType() == CbmTrdDigi::eCbmTrdAsicType::kSPADIC) { + vec->push_back(std::make_tuple(trdDigi->GetTime(), trdDigi->GetCharge(), trdDigi->GetAddressModule())); + return 1; } - catch (const std::invalid_argument& ia) { - if (fuNbTrdCondWarn < 10) { - LOG(warning) << "Trd condition for index " << uCond << " name " << det->vName[uCond] - << " not implemented. Skipped"; - ++fuNbTrdCondWarn; - } - return false; + + int toff(0); + uint n(0); + double t, r = trdDigi->GetCharge(t, toff); + if (t > 0) { + vec->push_back(std::make_tuple(dt, t, trdDigi->GetAddressModule())); + n++; } - return true; -} -bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmTofDigi* /*digi*/) -{ - return true; + if (r > 0) { + vec->push_back(std::make_tuple(dt + toff, r, trdDigi->GetAddressModule())); + n++; + } + return n; } -bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmRichDigi* /*digi*/) + +// ---- GetViewId ------------------------------------------------------ +template<class Digi> +int CbmMcbmCheckTimingAlgo::GetViewId(CheckTimingDetector, std::tuple<double, double, uint>) { - return true; + return 0; } -bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmPsdDigi* /*digi*/) +template<> +int CbmMcbmCheckTimingAlgo::GetViewId<CbmTrdDigi>(CheckTimingDetector det, std::tuple<double, double, uint> info) { - return true; + uint moduleId = std::get<2>(info); + uint iview = 0; + for (auto view : det.vName) { + if (view.compare(std::to_string(moduleId)) == 0) return iview; + iview++; + } + LOG(warning) << "Trd condition not implemented. Skip."; + return -1; } // ---- Finish -------------------------------------------------------- @@ -505,17 +517,17 @@ void CbmMcbmCheckTimingAlgo::WriteHistos() LOG(debug) << "Saving histos for " << uDet.sName; outfile->mkdir(uDet.sName.data()); outfile->cd(uDet.sName.data()); - for (uint id(0); id < uDet.uNdiv; id++) + for (uint id(0); id < uDet.uNviews; id++) fvhDetSelfDiff[uDet.detId][id]->Write(); - for (uint id(0); id < uDet.uNdiv; id++) + for (uint id(0); id < uDet.uNviews; id++) fvhDetToRefDiffRefCharge[uDet.detId][id]->Write(); - for (uint id(0); id < uDet.uNdiv; id++) + for (uint id(0); id < uDet.uNviews; id++) fvhDetToRefDiffDetCharge[uDet.detId][id]->Write(); - for (uint id(0); id < uDet.uNdiv; id++) + for (uint id(0); id < uDet.uNviews; id++) fvhDetToRefDiffEvo[uDet.detId][id]->Write(); - for (uint id(0); id < uDet.uNdiv; id++) + for (uint id(0); id < uDet.uNviews; id++) fvhDetToRefDiffEvoLong[uDet.detId][id]->Write(); - for (uint id(0); id < uDet.uNdiv; id++) { + for (uint id(0); id < uDet.uNviews; id++) { LOG(debug) << "WriteHistos, Det, entries = " << uDet.sName << " " << fvhDetToRefDiff[uDet.detId][id]->GetEntries(); LOG(info) << "Saved histos for " << uDet.sName << "(" << uDet.vName[id] << ")"; @@ -622,7 +634,7 @@ void CbmMcbmCheckTimingAlgo::WriteHistos() /// Register the histo for reference detector digi to digi time difference outfile->mkdir(fRefDet.sName.data()); outfile->cd(fRefDet.sName.data()); - for (uint id(0); id < fRefDet.uNdiv; id++) + for (uint id(0); id < fRefDet.uNviews; id++) fvhDetSelfDiff[fRefDet.detId][id]->Write(); outfile->cd(); @@ -655,7 +667,7 @@ void CbmMcbmCheckTimingAlgo::SetDetectorDifferential(ECbmModuleId detIn, std::ve if ((*det).detId != detIn) continue; found = true; (*det).vName = vName; - (*det).uNdiv = vName.size(); + (*det).uNviews = vName.size(); break; } @@ -671,15 +683,14 @@ void CbmMcbmCheckTimingAlgo::AddCheckDetector(ECbmModuleId detIn, std::string sN std::vector<CheckTimingDetector>::iterator det; for (det = fvDets.begin(); det != fvDets.end(); ++det) { if ((*det).detId == detIn) { - LOG(warning) << "CbmMcbmCheckTimingAlgo::AddCheckDetector => Detector already in " - "list, this call will only update the parameters!"; - (*det).dTimeRangeBeg = dTimeRangeBegIn; (*det).dTimeRangeEnd = dTimeRangeEndIn; (*det).uRangeNbBins = uRangeNbBinsIn; (*det).uChargeCutMin = uChargeCutMinIn; (*det).uChargeCutMax = uChargeCutMaxIn; - + LOG(info) << "CbmMcbmCheckTimingAlgo::AddCheckDetector => Detector " << (*det).sName << " range " + << (*det).dTimeRangeBeg << " :: " << (*det).dTimeRangeEnd << " bins " << (*det).uRangeNbBins + << " charge " << (*det).uChargeCutMin << "::" << (*det).uChargeCutMax; break; } // if( (*det).detId == detIn ) } // for( det = fvDets.begin(); det != fvDets.end(); ++det ) diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h index b4f4bc6f4cb2e0fedb2da717c632122d0d6283f2..65fb0ad61f8907c420d356a19f40eb8d4307c0b2 100644 --- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h +++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h @@ -42,8 +42,8 @@ public: 0; /// Charge cut used for example to reject/select pulser, no effect if equal, select if min < max, reject if max < min UInt_t uChargeCutMax = 0; /// Charge cut used for example to reject/select pulser, no effect if equal, select if min < max, reject if max < min - UInt_t uNdiv = 1; /// No of subdivisions for ech detector - std::vector<std::string> vName = {""}; + UInt_t uNviews = 1; /// No of views for each detector + std::vector<std::string> vName = {""}; /// view string definitions; to be defined by detectors /// Book-keeping variables Double_t dPrevTime = 0.; Int_t iPrevRefFirstDigi = 0; @@ -107,14 +107,22 @@ private: void CheckInterSystemOffset(); template<class Digi> void FillTimeOffsetHistos(const Double_t dRefTime, const Double_t dRefCharge, UInt_t uDetIdx); - bool CheckCondition(UInt_t uDet, UInt_t uCond, const void* digi); - bool CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmStsDigi* digi); - bool CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmMuchDigi* digi); - bool CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmTrdDigi* digi); - bool CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmTofDigi* digi); - bool CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmRichDigi* digi); - bool CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmPsdDigi* digi); - + /** @brief Retrieve digi (time,charge,addres) info. Use template specialization for special cases (e.g. T0, TRD, PSD) + * @param iDigi digi index in the digi array + * @param vec on return contains the signal(s), time(s) and address pairs. Should be allocated by the user + * @param detId if needed spec + * @return size of vector + */ + template<class Digi> + uint GetDigiInfo(UInt_t iDigi, std::vector<std::tuple<double, double, uint>>* vec, + ECbmModuleId detId = ECbmModuleId::kNotExist); + /** @brief Retrieve the detector view corresponding to the digi data (@see CheckTimingDetector::vName) + * @param det detector definitions + * @param info tuple of digi info (time, charge address) + * @return the view index for the curent data or -1 if there is none + */ + template<class Digi> + int GetViewId(CheckTimingDetector det, std::tuple<double, double, uint> info); /** Input array from previous already existing data level **/ CbmDigiManager* fDigiMan = nullptr; //! @@ -125,11 +133,10 @@ private: /** @brief Pointer to the Timeslice start time used to write it to the output tree @remark since we hand this to the FairRootManager it also wants to delete it and we do not have to take care of deletion **/ - CbmTsEventHeader* fCbmTsEventHeader = nullptr; + const CbmTsEventHeader* fCbmTsEventHeader = nullptr; // UInt_t fuNbTs = 0; - UInt_t fuNbTrdCondWarn = 0; CheckTimingDetector fRefDet {CheckTimingDetector(ECbmModuleId::kT0, "T0")}; std::vector<CheckTimingDetector> fvDets {