From ec8c80045c954a763fccda6809413ed6ef983c5c Mon Sep 17 00:00:00 2001 From: Alexandru Bercuci <abercuci@niham.nipne.ro> Date: Tue, 15 Mar 2022 23:07:40 +0000 Subject: [PATCH] extend use of check time difference algo for the case of user selections --- .../mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx | 427 +++++++++++------- fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h | 30 +- .../mcbm2018/tasks/CbmMcbmCheckTimingTask.cxx | 5 + fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.h | 1 + macro/beamtime/mcbm2022/check_timing_any.C | 143 ++++++ 5 files changed, 429 insertions(+), 177 deletions(-) create mode 100644 macro/beamtime/mcbm2022/check_timing_any.C diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx index ee6bacc49c..d878ead93d 100644 --- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx +++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.cxx @@ -97,9 +97,6 @@ void CbmMcbmCheckTimingAlgo::CheckDataPresence(CheckTimingDetector detToCheck) void CbmMcbmCheckTimingAlgo::CreateHistos() { - /// FIXME: Disable clang formatting for histograms declaration for now - /* clang-format off */ - /// Logarithmic bining for self time comparison uint32_t iNbBinsLog = 0; /// Parameters are NbDecadesLog, NbStepsDecade, NbSubStepsInStep @@ -108,51 +105,51 @@ void CbmMcbmCheckTimingAlgo::CreateHistos() for( std::vector< CheckTimingDetector >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) { - fvhDetSelfDiff.push_back( new TH1D( Form( "h%sSelfDiff", (*det).sName.data() ), - Form( "time difference between consecutivs %s Digis;time diff [ns];Counts", - (*det).sName.data() ), - iNbBinsLog, dBinsLog ) - ); - - fvhDetToRefDiff.push_back( new TH1D( Form( "h%s%sDiff", (*det).sName.data(), fRefDet.sName.data() ), - Form( "%s - %s time difference;time diff [ns];Counts", - (*det).sName.data(), fRefDet.sName.data() ), - (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd ) - ); - - fvhDetToRefDiffRefCharge.push_back( new TH2F( Form( "h%s%sDiffRefCharge", (*det).sName.data(), fRefDet.sName.data() ), - Form( "%s - %s;time diff [ns]; %s Charge [a.u]; Counts", - (*det).sName.data(), fRefDet.sName.data(), fRefDet.sName.data() ), - (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd, - 256, 0, 256 ) - ); - fvhDetToRefDiffDetCharge.push_back( new TH2F( Form( "h%s%sDiffDetCharge", (*det).sName.data(), fRefDet.sName.data() ), - Form( "%s - %s;time diff [ns]; %s Charge [a.u]; Counts", - (*det).sName.data(), fRefDet.sName.data(), (*det).sName.data() ), - (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd, - 256, 0, 256 ) - ); - fvhDetToRefDiffEvo.push_back( new TH2F( Form( "h%s%sDiffEvo", (*det).sName.data(), fRefDet.sName.data() ), - Form( "%s - %s;TS; time diff [ns];Counts", - (*det).sName.data(), fRefDet.sName.data() ), - 1000, 0, 10000, - (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd ) - ); - fvhDetToRefDiffEvoLong.push_back( new TH2F( Form( "h%s%sDiffEvoLong", (*det).sName.data(), fRefDet.sName.data() ), - Form( "%s - %s;TS; time diff [ns];Counts", - (*det).sName.data(), fRefDet.sName.data() ), - 1800, 0, 18000, - (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd ) - ); + for (uint i(0); i < (*det).uNdiv; 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(), + (*det).vName[i].data()), + iNbBinsLog, dBinsLog)); + + fvhDetToRefDiff[(*det).detId].push_back( + new TH1D(Form("h%s%sDiff%d", (*det).sName.data(), fRefDet.sName.data(), i), + Form("%s(%s) - %s time difference;time diff [ns];Counts", (*det).sName.data(), (*det).vName[i].data(), + fRefDet.sName.data()), + (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd)); + + fvhDetToRefDiffRefCharge[(*det).detId].push_back( + new TH2F(Form("h%s%sDiffRefCharge%d", (*det).sName.data(), fRefDet.sName.data(), i), + Form("%s(%s) - %s;time diff [ns]; %s Charge [a.u]; Counts", (*det).sName.data(), + (*det).vName[i].data(), fRefDet.sName.data(), fRefDet.sName.data()), + (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd, 256, 0, 256)); + + fvhDetToRefDiffDetCharge[(*det).detId].push_back( + new TH2F(Form("h%s%sDiffDetCharge%d", (*det).sName.data(), fRefDet.sName.data(), i), + Form("%s(%s) - %s;time diff [ns]; %s(%s) Charge [a.u]; Counts", (*det).sName.data(), + (*det).vName[i].data(), fRefDet.sName.data(), (*det).sName.data(), (*det).vName[i].data()), + (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd, 256, 0, 256)); + + fvhDetToRefDiffEvo[(*det).detId].push_back( + new TH2F(Form("h%s%sDiffEvo%d", (*det).sName.data(), fRefDet.sName.data(), i), + Form("%s(%s) - %s;TS; time diff [ns];Counts", (*det).sName.data(), (*det).vName[i].data(), + fRefDet.sName.data()), + 1000, 0, 10000, (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd)); + + fvhDetToRefDiffEvoLong[(*det).detId].push_back( + new TH2F(Form("h%s%sDiffEvoLong%d", (*det).sName.data(), fRefDet.sName.data(), i), + Form("%s(%s) - %s;TS; time diff [ns];Counts", (*det).sName.data(), (*det).vName[i].data(), + fRefDet.sName.data()), + 1800, 0, 18000, (*det).uRangeNbBins, (*det).dTimeRangeBeg, (*det).dTimeRangeEnd)); + } LOG( info ) << "Created histos for " << (*det).sName; } // for( std::vector< CheckTimingDetector >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) /// Add reference detector digi to digi time difference histo at end of vector - fvhDetSelfDiff.push_back( new TH1D( Form( "h%sSelfDiff", fRefDet.sName.data() ), - Form( "time difference between consecutivs %s Digis;time diff [ns];Counts", - fRefDet.sName.data() ), - iNbBinsLog, dBinsLog ) - ); + fvhDetSelfDiff[fRefDet.detId].push_back( + new TH1D(Form("h%sSelfDiff", fRefDet.sName.data()), + Form("time difference between consecutivs %s Digis;time diff [ns];Counts", fRefDet.sName.data()), + iNbBinsLog, dBinsLog)); /// Register the histos in the HTTP server FairRunOnline* run = FairRunOnline::Instance(); @@ -162,23 +159,19 @@ void CbmMcbmCheckTimingAlgo::CreateHistos() if( nullptr != server ) { /// Register histos for all checked detectors - for( UInt_t uDetIdx = 0; uDetIdx < fvDets.size(); ++uDetIdx ) - { - server->Register("/CheckTiming/SelfDiff", fvhDetSelfDiff[ uDetIdx ] ); - server->Register("/CheckTiming/RefDiff", fvhDetToRefDiff[ uDetIdx ] ); - server->Register("/CheckTiming/DiffCharge", fvhDetToRefDiffRefCharge[ uDetIdx ] ); - server->Register("/CheckTiming/DiffCharge", fvhDetToRefDiffDetCharge[ uDetIdx ] ); - server->Register("/CheckTiming/DiffEvo", fvhDetToRefDiffEvo[ uDetIdx ] ); - server->Register("/CheckTiming/DiffEvo", fvhDetToRefDiffEvoLong[ uDetIdx ] ); - } // for( std::vector< CheckTimingDetector >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) + for (auto uDetIdx : fvDets) { + server->Register("/CheckTiming/SelfDiff", fvhDetSelfDiff[uDetIdx.detId][0]); + server->Register("/CheckTiming/RefDiff", fvhDetToRefDiff[uDetIdx.detId][0]); + server->Register("/CheckTiming/DiffCharge", fvhDetToRefDiffRefCharge[uDetIdx.detId][0]); + server->Register("/CheckTiming/DiffCharge", fvhDetToRefDiffDetCharge[uDetIdx.detId][0]); + server->Register("/CheckTiming/DiffEvo", fvhDetToRefDiffEvo[uDetIdx.detId][0]); + server->Register("/CheckTiming/DiffEvo", fvhDetToRefDiffEvoLong[uDetIdx.detId][0]); + } // for( std::vector< CheckTimingDetector >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) /// Register the histo for reference detector digi to digi time difference - server->Register("/CheckTiming/SelfDiff", fvhDetSelfDiff[ fvDets.size() ] ); + server->Register("/CheckTiming/SelfDiff", fvhDetSelfDiff[fRefDet.detId][0]); } // if( nullptr != server ) } // if( run ) - - /// FIXME: Re-enable clang formatting after histograms declaration - /* clang-format on */ } // ---- ReInit ------------------------------------------------------- Bool_t CbmMcbmCheckTimingAlgo::ReInit() { return kTRUE; } @@ -272,7 +265,7 @@ void CbmMcbmCheckTimingAlgo::CheckInterSystemOffset() } // else of if( ECbmModuleId::kT0 == fRefDet.detId ) /// Fill self time difference histo - fvhDetSelfDiff[fvDets.size()]->Fill(dRefTime - fRefDet.dPrevTime); + (fvhDetSelfDiff[fRefDet.detId])[0]->Fill(dRefTime - fRefDet.dPrevTime); fRefDet.dPrevTime = dRefTime; /// Charge cut if defined! @@ -342,7 +335,8 @@ template<class Digi> void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const Double_t dRefCharge, UInt_t uDetIdx) { UInt_t uNbDigis = 0; - switch (fvDets[uDetIdx].detId) { + ECbmModuleId edetId = fvDets[uDetIdx].detId; + switch (edetId) { case ECbmModuleId::kNotExist: { LOG(fatal) << "CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos => Unknow " "detector enum! " @@ -354,22 +348,25 @@ void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const break; } // case ECbmModuleId::kT0 default: { - uNbDigis = fDigiMan->GetNofDigis(fvDets[uDetIdx].detId); + uNbDigis = fDigiMan->GetNofDigis(edetId); break; } // default: } // switch( fRefDet.detId ) UInt_t uFirstDigiInWin = fvDets[uDetIdx].iPrevRefFirstDigi; + std::vector<Double_t> vSelDiff; for (UInt_t uDigiIdx = fvDets[uDetIdx].iPrevRefFirstDigi; uDigiIdx < uNbDigis; ++uDigiIdx) { Double_t dTime = 0; Double_t dCharge = 0; UInt_t uAddress = 0; - if (ECbmModuleId::kT0 == fvDets[uDetIdx].detId) { + 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(); @@ -378,7 +375,7 @@ void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const /// Fill self correlation histo while avoiding double counting due to /// the "smart looping" if (fvDets[uDetIdx].dPrevTime <= dTime) { - fvhDetSelfDiff[uDetIdx]->Fill(dTime - fvDets[uDetIdx].dPrevTime); + vSelDiff.push_back(dTime - fvDets[uDetIdx].dPrevTime); fvDets[uDetIdx].dPrevTime = dTime; } // if( fvDets[ uDetIdx ].dPrevTime < dTime ) @@ -411,17 +408,26 @@ void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const } // else of if( fvDets[ uDetIdx ].uChargeCutMin < fvDets[ uDetIdx ].uChargeCutMax ) } // if( fvDets[ uDetIdx ].uChargeCutMin != fvDets[ uDetIdx ].uChargeCutMax ) - /// Fill histos - fvhDetToRefDiff[uDetIdx]->Fill(dDiffTime); - fvhDetToRefDiffRefCharge[uDetIdx]->Fill(dDiffTime, dRefCharge); - fvhDetToRefDiffDetCharge[uDetIdx]->Fill(dDiffTime, dCharge); - if (nullptr == fCbmTsEventHeader) { - fvhDetToRefDiffEvo[uDetIdx]->Fill(fuNbTs, dDiffTime); - fvhDetToRefDiffEvoLong[uDetIdx]->Fill(fuNbTs, dDiffTime); - } - else { - fvhDetToRefDiffEvo[uDetIdx]->Fill(fCbmTsEventHeader->GetTsIndex(), dDiffTime); - fvhDetToRefDiffEvoLong[uDetIdx]->Fill(fCbmTsEventHeader->GetTsIndex(), dDiffTime); + for (uint uid(0); uid < fvDets[uDetIdx].uNdiv; uid++) { + if (digi && !CheckCondition(uDetIdx, uid, digi)) continue; + + /// Fill histos + for (auto dt : vSelDiff) + fvhDetSelfDiff[edetId][uid]->Fill(dt); + vSelDiff.clear(); + + fvhDetToRefDiff[edetId][uid]->Fill(dDiffTime); + fvhDetToRefDiffRefCharge[edetId][uid]->Fill(dDiffTime, dRefCharge); + fvhDetToRefDiffDetCharge[edetId][uid]->Fill(dDiffTime, dCharge); + if (nullptr == fCbmTsEventHeader) { + fvhDetToRefDiffEvo[edetId][uid]->Fill(fuNbTs, dDiffTime); + fvhDetToRefDiffEvoLong[edetId][uid]->Fill(fuNbTs, dDiffTime); + } + else { + fvhDetToRefDiffEvo[edetId][uid]->Fill(fCbmTsEventHeader->GetTsIndex(), dDiffTime); + fvhDetToRefDiffEvoLong[edetId][uid]->Fill(fCbmTsEventHeader->GetTsIndex(), dDiffTime); + } + break; } } // for( UInt_t uDigiIdx = fvDets[ uDetIdx ].iPrevRefFirstDigi; uDigiIdx < uNbDigis; ++uDigiIdx ) @@ -429,6 +435,55 @@ void CbmMcbmCheckTimingAlgo::FillTimeOffsetHistos(const Double_t dRefTime, const 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*/) +{ + return true; +} +bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmMuchDigi* /*digi*/) +{ + return true; +} +bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* det, UInt_t uCond, const CbmTrdDigi* digi) +{ + try { + int moduleId = std::stoi(det->vName[uCond]); + if (digi->GetAddressModule() != moduleId) return false; + } + catch (const std::invalid_argument& ia) { + LOG(warning) << "Trd condition " << det->vName[uCond] << " not implemented. Skipped"; + return false; + } + return true; +} +bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmTofDigi* /*digi*/) +{ + return true; +} +bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmRichDigi* /*digi*/) +{ + return true; +} +bool CbmMcbmCheckTimingAlgo::CheckCondition(CheckTimingDetector* /*det*/, UInt_t /*uCond*/, const CbmPsdDigi* /*digi*/) +{ + return true; +} + // ---- Finish -------------------------------------------------------- void CbmMcbmCheckTimingAlgo::Finish() { LOG(info) << Form("Checked %6d Timeslices", fuNbTs); } @@ -439,114 +494,130 @@ void CbmMcbmCheckTimingAlgo::WriteHistos() TFile* outfile = TFile::Open(fOutFileName, "RECREATE"); - for (UInt_t uDetIdx = 0; uDetIdx < fvDets.size(); ++uDetIdx) { - LOG(debug) << "Saving histos for " << fvDets[uDetIdx].sName; - fvhDetSelfDiff[uDetIdx]->Write(); - fvhDetToRefDiffRefCharge[uDetIdx]->Write(); - fvhDetToRefDiffDetCharge[uDetIdx]->Write(); - fvhDetToRefDiffEvo[uDetIdx]->Write(); - fvhDetToRefDiffEvoLong[uDetIdx]->Write(); - LOG(debug) << "WriteHistos, uDetIdx, Det, entries = " << uDetIdx << " " << fvDets[uDetIdx].sName << " " - << fvhDetToRefDiff[uDetIdx]->GetEntries(); - LOG(info) << "Saved histos for " << fvDets[uDetIdx].sName; - DetPeakPosSingle = fvhDetToRefDiff[uDetIdx]->GetMaximumBin() * fvhDetToRefDiff[uDetIdx]->GetBinWidth(1) - + fvhDetToRefDiff[uDetIdx]->GetXaxis()->GetXmin(); - DetAverageSingle = (fvhDetToRefDiff[uDetIdx]->Integral()) / (fvhDetToRefDiff[uDetIdx]->GetNbinsX()); - switch (fvDets[uDetIdx].detId) { - case ECbmModuleId::kSts: { - if (DetAverageSingle > 0) { - TF1* gs_sts = new TF1("gs_sts", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fStsPeakWidthNs, - DetPeakPosSingle + 2 * fStsPeakWidthNs); - gs_sts->SetParameters(DetAverageSingle, DetPeakPosSingle, fStsPeakWidthNs, DetAverageSingle); - fvhDetToRefDiff[uDetIdx]->Fit("gs_sts", "R"); - TF1* fitresult_sts = fvhDetToRefDiff[uDetIdx]->GetFunction("gs_sts"); - LOG(debug) << fvDets[uDetIdx].sName << " parameters from Gauss fit = " << fitresult_sts->GetParameter(0) - << ", " << fitresult_sts->GetParameter(1) << ", " << fitresult_sts->GetParameter(2); + for (auto uDet : fvDets) { + 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++) + fvhDetSelfDiff[uDet.detId][id]->Write(); + for (uint id(0); id < uDet.uNdiv; id++) + fvhDetToRefDiffRefCharge[uDet.detId][id]->Write(); + for (uint id(0); id < uDet.uNdiv; id++) + fvhDetToRefDiffDetCharge[uDet.detId][id]->Write(); + for (uint id(0); id < uDet.uNdiv; id++) + fvhDetToRefDiffEvo[uDet.detId][id]->Write(); + for (uint id(0); id < uDet.uNdiv; id++) + fvhDetToRefDiffEvoLong[uDet.detId][id]->Write(); + for (uint id(0); id < uDet.uNdiv; id++) { + LOG(debug) << "WriteHistos, Det, entries = " << uDet.sName << " " + << fvhDetToRefDiff[uDet.detId][id]->GetEntries(); + LOG(info) << "Saved histos for " << uDet.sName << "(" << uDet.vName[id] << ")"; + DetPeakPosSingle = + fvhDetToRefDiff[uDet.detId][id]->GetMaximumBin() * fvhDetToRefDiff[uDet.detId][id]->GetBinWidth(1) + + fvhDetToRefDiff[uDet.detId][id]->GetXaxis()->GetXmin(); + DetAverageSingle = (fvhDetToRefDiff[uDet.detId][id]->Integral()) / (fvhDetToRefDiff[uDet.detId][id]->GetNbinsX()); + + switch (uDet.detId) { + case ECbmModuleId::kSts: { + if (DetAverageSingle > 0) { + TF1* gs_sts = new TF1("gs_sts", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fStsPeakWidthNs, + DetPeakPosSingle + 2 * fStsPeakWidthNs); + gs_sts->SetParameters(DetAverageSingle, DetPeakPosSingle, fStsPeakWidthNs, DetAverageSingle); + fvhDetToRefDiff[uDet.detId][id]->Fit("gs_sts", "R"); + TF1* fitresult_sts = fvhDetToRefDiff[uDet.detId][id]->GetFunction("gs_sts"); + LOG(debug) << uDet.sName << " parameters from Gauss fit = " << fitresult_sts->GetParameter(0) << ", " + << fitresult_sts->GetParameter(1) << ", " << fitresult_sts->GetParameter(2); + } + break; } - break; - } - case ECbmModuleId::kMuch: { - if (DetAverageSingle > 0) { - TF1* gs_much = new TF1("gs_much", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fMuchPeakWidthNs, - DetPeakPosSingle + 2 * fMuchPeakWidthNs); - gs_much->SetParameters(DetAverageSingle, DetPeakPosSingle, fMuchPeakWidthNs, DetAverageSingle); - fvhDetToRefDiff[uDetIdx]->Fit("gs_much", "R"); - TF1* fitresult_much = fvhDetToRefDiff[uDetIdx]->GetFunction("gs_much"); - LOG(debug) << fvDets[uDetIdx].sName << " parameters from Gauss fit = " << fitresult_much->GetParameter(0) - << ", " << fitresult_much->GetParameter(1) << ", " << fitresult_much->GetParameter(2); + case ECbmModuleId::kMuch: { + if (DetAverageSingle > 0) { + TF1* gs_much = new TF1("gs_much", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fMuchPeakWidthNs, + DetPeakPosSingle + 2 * fMuchPeakWidthNs); + gs_much->SetParameters(DetAverageSingle, DetPeakPosSingle, fMuchPeakWidthNs, DetAverageSingle); + fvhDetToRefDiff[uDet.detId][id]->Fit("gs_much", "R"); + TF1* fitresult_much = fvhDetToRefDiff[uDet.detId][id]->GetFunction("gs_much"); + LOG(debug) << uDet.sName << " parameters from Gauss fit = " << fitresult_much->GetParameter(0) << ", " + << fitresult_much->GetParameter(1) << ", " << fitresult_much->GetParameter(2); + } + break; } - break; - } - case ECbmModuleId::kTrd: { - if (DetAverageSingle > 0) { - TF1* gs_trd = new TF1("gs_trd", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fTrdPeakWidthNs, - DetPeakPosSingle + 2 * fTrdPeakWidthNs); - gs_trd->SetParameters(0.7 * DetAverageSingle, DetPeakPosSingle, fTrdPeakWidthNs, DetAverageSingle); - fvhDetToRefDiff[uDetIdx]->Fit("gs_trd", "R"); - TF1* fitresult_trd = fvhDetToRefDiff[uDetIdx]->GetFunction("gs_trd"); - LOG(debug) << fvDets[uDetIdx].sName << " parameters from Gauss fit = " << fitresult_trd->GetParameter(0) - << ", " << fitresult_trd->GetParameter(1) << ", " << fitresult_trd->GetParameter(2); + case ECbmModuleId::kTrd: { + if (DetAverageSingle > 0) { + TF1* gs_trd = new TF1("gs_trd", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fTrdPeakWidthNs, + DetPeakPosSingle + 2 * fTrdPeakWidthNs); + gs_trd->SetParameters(0.7 * DetAverageSingle, DetPeakPosSingle, fTrdPeakWidthNs, DetAverageSingle); + fvhDetToRefDiff[uDet.detId][id]->Fit("gs_trd", "R"); + TF1* fitresult_trd = fvhDetToRefDiff[uDet.detId][id]->GetFunction("gs_trd"); + LOG(debug) << uDet.sName << " parameters from Gauss fit = " << fitresult_trd->GetParameter(0) << ", " + << fitresult_trd->GetParameter(1) << ", " << fitresult_trd->GetParameter(2); + } + break; } - break; - } - case ECbmModuleId::kT0: { - if (DetAverageSingle > 0) { - TF1* gs_tof = new TF1("gs_tof", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fTofPeakWidthNs, - DetPeakPosSingle + 2 * fTofPeakWidthNs); - gs_tof->SetParameters(DetAverageSingle, DetPeakPosSingle, fTofPeakWidthNs, DetAverageSingle); - fvhDetToRefDiff[uDetIdx]->Fit("gs_tof", "R"); - TF1* fitresult_tof = fvhDetToRefDiff[uDetIdx]->GetFunction("gs_tof"); - LOG(debug) << fvDets[uDetIdx].sName << " parameters from Gauss fit = " << fitresult_tof->GetParameter(0) - << ", " << fitresult_tof->GetParameter(1) << ", " << fitresult_tof->GetParameter(2); + case ECbmModuleId::kT0: { + if (DetAverageSingle > 0) { + TF1* gs_tof = new TF1("gs_tof", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fTofPeakWidthNs, + DetPeakPosSingle + 2 * fTofPeakWidthNs); + gs_tof->SetParameters(DetAverageSingle, DetPeakPosSingle, fTofPeakWidthNs, DetAverageSingle); + fvhDetToRefDiff[uDet.detId][id]->Fit("gs_tof", "R"); + TF1* fitresult_tof = fvhDetToRefDiff[uDet.detId][id]->GetFunction("gs_tof"); + LOG(debug) << uDet.sName << " parameters from Gauss fit = " << fitresult_tof->GetParameter(0) << ", " + << fitresult_tof->GetParameter(1) << ", " << fitresult_tof->GetParameter(2); + } + break; } - break; - } - case ECbmModuleId::kTof: { - if (DetAverageSingle > 0) { - TF1* gs_tof = new TF1("gs_tof", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fTofPeakWidthNs, - DetPeakPosSingle + 2 * fTofPeakWidthNs); - gs_tof->SetParameters(DetAverageSingle, DetPeakPosSingle, fTofPeakWidthNs, DetAverageSingle); - fvhDetToRefDiff[uDetIdx]->Fit("gs_tof", "R"); - TF1* fitresult_tof = fvhDetToRefDiff[uDetIdx]->GetFunction("gs_tof"); - LOG(debug) << fvDets[uDetIdx].sName << " parameters from Gauss fit = " << fitresult_tof->GetParameter(0) - << ", " << fitresult_tof->GetParameter(1) << ", " << fitresult_tof->GetParameter(2); + case ECbmModuleId::kTof: { + if (DetAverageSingle > 0) { + TF1* gs_tof = new TF1("gs_tof", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fTofPeakWidthNs, + DetPeakPosSingle + 2 * fTofPeakWidthNs); + gs_tof->SetParameters(DetAverageSingle, DetPeakPosSingle, fTofPeakWidthNs, DetAverageSingle); + fvhDetToRefDiff[uDet.detId][id]->Fit("gs_tof", "R"); + TF1* fitresult_tof = fvhDetToRefDiff[uDet.detId][id]->GetFunction("gs_tof"); + LOG(debug) << uDet.sName << " parameters from Gauss fit = " << fitresult_tof->GetParameter(0) << ", " + << fitresult_tof->GetParameter(1) << ", " << fitresult_tof->GetParameter(2); + } + break; } - break; - } - case ECbmModuleId::kRich: { - if (DetAverageSingle > 0) { - TF1* gs_rich = new TF1("gs_rich", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fRichPeakWidthNs, - DetPeakPosSingle + 2 * fRichPeakWidthNs); - gs_rich->SetParameters(0.5 * DetAverageSingle, DetPeakPosSingle, fRichPeakWidthNs, DetAverageSingle); - fvhDetToRefDiff[uDetIdx]->Fit("gs_rich", "R"); - TF1* fitresult_rich = fvhDetToRefDiff[uDetIdx]->GetFunction("gs_rich"); - LOG(debug) << fvDets[uDetIdx].sName << " parameters from Gauss fit = " << fitresult_rich->GetParameter(0) - << ", " << fitresult_rich->GetParameter(1) << ", " << fitresult_rich->GetParameter(2); + case ECbmModuleId::kRich: { + if (DetAverageSingle > 0) { + TF1* gs_rich = new TF1("gs_rich", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fRichPeakWidthNs, + DetPeakPosSingle + 2 * fRichPeakWidthNs); + gs_rich->SetParameters(0.5 * DetAverageSingle, DetPeakPosSingle, fRichPeakWidthNs, DetAverageSingle); + fvhDetToRefDiff[uDet.detId][id]->Fit("gs_rich", "R"); + TF1* fitresult_rich = fvhDetToRefDiff[uDet.detId][id]->GetFunction("gs_rich"); + LOG(debug) << uDet.sName << " parameters from Gauss fit = " << fitresult_rich->GetParameter(0) << ", " + << fitresult_rich->GetParameter(1) << ", " << fitresult_rich->GetParameter(2); + } + break; } - break; - } - case ECbmModuleId::kPsd: { - if (DetAverageSingle > 0) { - TF1* gs_psd = new TF1("gs_psd", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fPsdPeakWidthNs, - DetPeakPosSingle + 2 * fPsdPeakWidthNs); - gs_psd->SetParameters(DetAverageSingle, DetPeakPosSingle, fPsdPeakWidthNs, DetAverageSingle); - fvhDetToRefDiff[uDetIdx]->Fit("gs_psd", "R"); - TF1* fitresult_psd = fvhDetToRefDiff[uDetIdx]->GetFunction("gs_psd"); - LOG(debug) << fvDets[uDetIdx].sName << " parameters from Gauss fit = " << fitresult_psd->GetParameter(0) - << ", " << fitresult_psd->GetParameter(1) << ", " << fitresult_psd->GetParameter(2); + case ECbmModuleId::kPsd: { + if (DetAverageSingle > 0) { + TF1* gs_psd = new TF1("gs_psd", "gaus(0)+pol0(3)", DetPeakPosSingle - 2 * fPsdPeakWidthNs, + DetPeakPosSingle + 2 * fPsdPeakWidthNs); + gs_psd->SetParameters(DetAverageSingle, DetPeakPosSingle, fPsdPeakWidthNs, DetAverageSingle); + fvhDetToRefDiff[uDet.detId][id]->Fit("gs_psd", "R"); + TF1* fitresult_psd = fvhDetToRefDiff[uDet.detId][id]->GetFunction("gs_psd"); + LOG(debug) << uDet.sName << " parameters from Gauss fit = " << fitresult_psd->GetParameter(0) << ", " + << fitresult_psd->GetParameter(1) << ", " << fitresult_psd->GetParameter(2); + } + break; + } + default: { + LOG(info) << "Detector ID for fitting is not valid."; + break; } - break; - } - default: { - LOG(info) << "Detector ID for fitting is not valid."; - break; } + fvhDetToRefDiff[uDet.detId][id]->Write(); //At the end in order to include fitting results in histos } - fvhDetToRefDiff[uDetIdx]->Write(); //At the end in order to include fitting results in histos + outfile->cd(); } // for( std::vector< CheckTimingDetector >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) /// Register the histo for reference detector digi to digi time difference - fvhDetSelfDiff[fvDets.size()]->Write(); + outfile->mkdir(fRefDet.sName.data()); + outfile->cd(fRefDet.sName.data()); + for (uint id(0); id < fRefDet.uNdiv; id++) + fvhDetSelfDiff[fRefDet.detId][id]->Write(); + outfile->cd(); outfile->Close(); delete outfile; @@ -568,6 +639,24 @@ void CbmMcbmCheckTimingAlgo::SetReferenceDetector(ECbmModuleId refDetIn, std::st fRefDet.uChargeCutMin = uChargeCutMinIn; fRefDet.uChargeCutMax = uChargeCutMaxIn; } + +void CbmMcbmCheckTimingAlgo::SetDetectorDifferential(ECbmModuleId detIn, std::vector<std::string> vName) +{ + bool found(false); + std::vector<CheckTimingDetector>::iterator det; + for (det = fvDets.begin(); det != fvDets.end(); ++det) { + if ((*det).detId != detIn) continue; + found = true; + (*det).vName = vName; + (*det).uNdiv = vName.size(); + break; + } + + if (!found) + LOG(warning) << "CbmMcbmCheckTimingAlgo::SetDetectorDifferential => Detector not in the " + "list, Nothing done at this point!"; +} + void CbmMcbmCheckTimingAlgo::AddCheckDetector(ECbmModuleId detIn, std::string sNameIn, Double_t dTimeRangeBegIn, Double_t dTimeRangeEndIn, UInt_t uRangeNbBinsIn, UInt_t uChargeCutMinIn, UInt_t uChargeCutMaxIn) diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h index 0ca30a3c6d..f1e484f999 100644 --- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h +++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingAlgo.h @@ -17,7 +17,12 @@ class TH1; class TH2; class CbmDigiManager; - +class CbmStsDigi; +class CbmMuchDigi; +class CbmTrdDigi; +class CbmTofDigi; +class CbmRichDigi; +class CbmPsdDigi; class CheckTimingDetector { public: CheckTimingDetector() { ; } @@ -37,7 +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 = {"none"}; /// Book-keeping variables Double_t dPrevTime = 0.; Int_t iPrevRefFirstDigi = 0; @@ -84,6 +90,7 @@ public: Double_t dTimeRangeEndIn = 1000.0, UInt_t uRangeNbBinsIn = 320, UInt_t uChargeCutMinIn = 0, UInt_t uChargeCutMaxIn = 0); void RemoveCheckDetector(ECbmModuleId detIn); + void SetDetectorDifferential(ECbmModuleId detIn, std::vector<std::string> vName); void SetTrdPeakWidthNs(Double_t val = 120.) { fTrdPeakWidthNs = val; } void SetStsPeakWidthNs(Double_t val = 30.) { fStsPeakWidthNs = val; } @@ -100,6 +107,13 @@ 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); /** Input array from previous already existing data level **/ @@ -123,12 +137,12 @@ private: CheckTimingDetector(ECbmModuleId::kRich, "Rich"), CheckTimingDetector(ECbmModuleId::kPsd, "Psd")}; /// vectors storing histograms for each detector under investigation - std::vector<TH1*> fvhDetSelfDiff = {}; - std::vector<TH1*> fvhDetToRefDiff = {}; - std::vector<TH2*> fvhDetToRefDiffRefCharge = {}; - std::vector<TH2*> fvhDetToRefDiffDetCharge = {}; - std::vector<TH2*> fvhDetToRefDiffEvo = {}; - std::vector<TH2*> fvhDetToRefDiffEvoLong = {}; + std::map<ECbmModuleId, std::vector<TH1*>> fvhDetSelfDiff = {}; + std::map<ECbmModuleId, std::vector<TH1*>> fvhDetToRefDiff = {}; + std::map<ECbmModuleId, std::vector<TH2*>> fvhDetToRefDiffRefCharge = {}; + std::map<ECbmModuleId, std::vector<TH2*>> fvhDetToRefDiffDetCharge = {}; + std::map<ECbmModuleId, std::vector<TH2*>> fvhDetToRefDiffEvo = {}; + std::map<ECbmModuleId, std::vector<TH2*>> fvhDetToRefDiffEvoLong = {}; /** Name of the histogram output file **/ TString fOutFileName = "data/HistosCheckTiming.root"; diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.cxx b/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.cxx index 52ef5cac8a..0831cbf6ae 100644 --- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.cxx +++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.cxx @@ -128,5 +128,10 @@ void CbmMcbmCheckTimingTask::AddCheckDetector(ECbmModuleId detIn, std::string sN void CbmMcbmCheckTimingTask::RemoveCheckDetector(ECbmModuleId detIn) { fpAlgo->RemoveCheckDetector(detIn); } +void CbmMcbmCheckTimingTask::SetDetectorDifferential(ECbmModuleId detIn, std::vector<std::string> vn) +{ + fpAlgo->SetDetectorDifferential(detIn, vn); +} + //---------------------------------------------------------------------- ClassImp(CbmMcbmCheckTimingTask) diff --git a/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.h b/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.h index 2ede5ce42a..b756ce0b50 100644 --- a/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.h +++ b/fles/mcbm2018/tasks/CbmMcbmCheckTimingTask.h @@ -58,6 +58,7 @@ public: Double_t dTimeRangeEndIn = 1000.0, UInt_t uRangeNbBinsIn = 320, UInt_t uChargeCutMinIn = 0, UInt_t uChargeCutMaxIn = 0); void RemoveCheckDetector(ECbmModuleId detIn); + void SetDetectorDifferential(ECbmModuleId detIn, std::vector<std::string> vName); private: void SaveHistos(); diff --git a/macro/beamtime/mcbm2022/check_timing_any.C b/macro/beamtime/mcbm2022/check_timing_any.C new file mode 100644 index 0000000000..7535424412 --- /dev/null +++ b/macro/beamtime/mcbm2022/check_timing_any.C @@ -0,0 +1,143 @@ +/* Copyright (C) 2020 Facility for Antiproton and Ion Research in Europe, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Pierre-Alain Loizeau [committer] */ + +void check_timing_any(TString fileName, UInt_t uRunId = 0, Int_t nEvents = 0, TString outDir = "data/") +{ + + // ======================================================================== + // Adjust this part according to your requirements + + // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) + Int_t iVerbose = 3; + + // MC file + + TString srcDir = gSystem->Getenv("VMCWORKDIR"); + + // ----- Timer -------------------------------------------------------- + TStopwatch timer; + timer.Start(); + // ------------------------------------------------------------------------ + + // ----- Analysis run -------------------------------------------------- + FairRunOnline* fRun = new FairRunOnline(); + fRun->ActivateHttpServer(100, 8080); // refresh each 100 events + fRun->SetSink(new FairRootFileSink("SinkFile.root")); + FairFileSource* inputSource = new FairFileSource(fileName); + fRun->SetSource(inputSource); + + // Define output file for FairMonitor histograms + // TString monitorFile{outFile}; + // monitorFile.ReplaceAll("qa","qa.monitor"); + FairMonitor::GetMonitor()->EnableMonitor(kFALSE); + // ------------------------------------------------------------------------ + + CbmMcbmCheckTimingTask* timeChecker = new CbmMcbmCheckTimingTask(); + /// Default is using T0 as reference + /// With Pulser rejection + /* + timeChecker->SetReferenceDetector( ECbmModuleId::kT0, "T0", + -1000., 1000., 320., + 182, 190 ); +*/ + /// With pulser selection + /* + timeChecker->SetReferenceDetector( ECbmModuleId::kT0, "T0", + -1000., 1000., 320., + 190, 182 ); +*/ + /// Here swapping with TOF + timeChecker->SetReferenceDetector(ECbmModuleId::kTof, "Tof", -3000, 3000, 3000); + timeChecker->RemoveCheckDetector(ECbmModuleId::kTof); + //timeChecker->AddCheckDetector(ECbmModuleId::kT0, "T0"); + + /// Here swapping with MUCH + /* + timeChecker->SetReferenceDetector(ECbmModuleId::kMuch, "Much"); + timeChecker->RemoveCheckDetector(ECbmModuleId::kMuch); + timeChecker->AddCheckDetector(ECbmModuleId::kT0, "T0"); + */ + + /// Remove detectors not present in 2021 + timeChecker->RemoveCheckDetector(ECbmModuleId::kT0); + timeChecker->RemoveCheckDetector(ECbmModuleId::kMuch); + //timeChecker->SetReferenceDetector(ECbmModuleId::kPsd, "Psd", -300000, 300000, 320 * 300); + timeChecker->RemoveCheckDetector(ECbmModuleId::kPsd); + + /// Add detectors with wider range + /* + timeChecker->RemoveCheckDetector(ECbmModuleId::kSts); + timeChecker->AddCheckDetector(ECbmModuleId::kSts, "Sts"); + timeChecker->RemoveCheckDetector(ECbmModuleId::kTrd); + timeChecker->AddCheckDetector(ECbmModuleId::kTrd, "Trd"); + timeChecker->RemoveCheckDetector(ECbmModuleId::kTof); + //timeChecker->AddCheckDetector(ECbmModuleId::kTof, "Tof", -150000, 150000, 320*150); + timeChecker->AddCheckDetector(ECbmModuleId::kTof, "Tof", -2000, 2000, 320 * 2); + */ + + /// Add extra differential analysis on specific detectors + timeChecker->AddCheckDetector(ECbmModuleId::kSts, "Sts", -200., 200., 4000); + timeChecker->AddCheckDetector(ECbmModuleId::kRich, "Rich", -100., 100., 4000); + timeChecker->AddCheckDetector(ECbmModuleId::kTrd, "Trd", -500.5, 500.5, 1001); + timeChecker->SetDetectorDifferential(ECbmModuleId::kTrd, {"5", "21", "37", "53"}); + + TString outFileName = fileName(0, fileName.Index(".digi.root")); + if (0 < uRunId) timeChecker->SetOutFilename(Form("%s/%s.tck.root", outDir.Data(), outFileName.Data())); + fRun->AddTask(timeChecker); + + // ----- Parameter database -------------------------------------------- + // FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); + // FairParRootFileIo* parIo1 = new FairParRootFileIo(); + // parIo1->open(parFile.Data(),"UPDATE"); + // rtdb->setFirstInput(parIo1); + // ------------------------------------------------------------------------ + + + // ----- Intialise and run -------------------------------------------- + fRun->Init(); + + // rtdb->setOutput(parIo1); + // rtdb->saveOutput(); + // rtdb->print(); + + cout << "Starting run" << endl; + if (0 == nEvents) { + fRun->Run(0, 0); // run until end of input file + } + else { + fRun->Run(0, nEvents); // process N Events + } + // ------------------------------------------------------------------------ + + + // ----- Finish ------------------------------------------------------- + timer.Stop(); + Double_t rtime = timer.RealTime(); + Double_t ctime = timer.CpuTime(); + cout << endl << endl; + cout << "Macro finished successfully." << endl; + cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl; + cout << endl; + // ------------------------------------------------------------------------ + + // Extract the maximal used memory an add is as Dart measurement + // This line is filtered by CTest and the value send to CDash + FairSystemInfo sysInfo; + Float_t maxMemory = sysInfo.GetMaxMemory(); + cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">"; + cout << maxMemory; + cout << "</DartMeasurement>" << endl; + + Float_t cpuUsage = ctime / rtime; + cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">"; + cout << cpuUsage; + cout << "</DartMeasurement>" << endl; + /* + FairMonitor* tempMon = FairMonitor::GetMonitor(); + tempMon->Print(); +*/ + // RemoveGeoManager(); + cout << " Test passed" << endl; + cout << " All ok " << endl; +} -- GitLab