diff --git a/fles/mcbm2018/CMakeLists.txt b/fles/mcbm2018/CMakeLists.txt index e44e2b2cc4523712171294ce8a5d36649bb82d3b..dd306e6d230d303f655a555b4b6efd4458b74e7b 100644 --- a/fles/mcbm2018/CMakeLists.txt +++ b/fles/mcbm2018/CMakeLists.txt @@ -122,6 +122,10 @@ Set(SRCS tasks/CbmMcbm2019CheckPulser.cxx tasks/CbmMcbm2019CheckDigisMuch.cxx tasks/CbmMcbm2019CheckDigisSts.cxx + tasks/CbmMcbm2019CheckDtInDet.cxx + tasks/CbmMcbm2019CheckTimingPairs.cxx + tasks/CbmMcbm2019TimeWinEventBuilderAlgo.cxx + tasks/CbmMcbm2019TimeWinEventBuilderTask.cxx ) If(_UINT8_T_EXIST) diff --git a/fles/mcbm2018/CbmFlibMcbm2018LinkDef.h b/fles/mcbm2018/CbmFlibMcbm2018LinkDef.h index f0c8469becde5a604925a6ae4518408ec3a33135..3276e140224d56d805bc570297fbcbbba5f7c055 100644 --- a/fles/mcbm2018/CbmFlibMcbm2018LinkDef.h +++ b/fles/mcbm2018/CbmFlibMcbm2018LinkDef.h @@ -68,5 +68,9 @@ #pragma link C++ class CbmMcbm2019CheckPulser; #pragma link C++ class CbmMcbm2019CheckDigisMuch; #pragma link C++ class CbmMcbm2019CheckDigisSts; +#pragma link C++ class CbmMcbm2019CheckDtInDet; +#pragma link C++ class CbmMcbm2019CheckTimingPairs+; +#pragma link C++ class CbmMcbm2019TimeWinEventBuilderAlgo+; +#pragma link C++ class CbmMcbm2019TimeWinEventBuilderTask+; #endif diff --git a/fles/mcbm2018/tasks/CbmMcbm2019CheckDtInDet.cxx b/fles/mcbm2018/tasks/CbmMcbm2019CheckDtInDet.cxx new file mode 100644 index 0000000000000000000000000000000000000000..bb3e4276611d2f068056c41cc42dfc68903f2a49 --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019CheckDtInDet.cxx @@ -0,0 +1,474 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#include "CbmMcbm2019CheckDtInDet.h" + +#include "CbmStsDigi.h" +#include "CbmMuchBeamTimeDigi.h" +#include "CbmRichDigi.h" +#include "CbmTrdDigi.h" +#include "CbmTofDigi.h" +#include "CbmPsdDigi.h" +#include "CbmDigiManager.h" +#include "CbmFlesHistosTools.h" + +#include "FairLogger.h" +#include "FairRootManager.h" +#include "FairRunOnline.h" + +#include "TClonesArray.h" +#include "TH1.h" +#include "TH2.h" +#include "TProfile.h" +#include "THttpServer.h" +#include <TFile.h> + + +#include <iomanip> +#include <iostream> +#include <type_traits> +using std::fixed; +using std::setprecision; + +// ---- Default constructor ------------------------------------------- +CbmMcbm2019CheckDtInDet::CbmMcbm2019CheckDtInDet() + : FairTask("CbmMcbm2019CheckDtInDet") +{ +} + +// ---- Destructor ---------------------------------------------------- +CbmMcbm2019CheckDtInDet::~CbmMcbm2019CheckDtInDet() +{ +} + +// ---- Initialisation ---------------------------------------------- +void CbmMcbm2019CheckDtInDet::SetParContainers() +{ + // Load all necessary parameter containers from the runtime data base + /* + FairRunAna* ana = FairRunAna::Instance(); + FairRuntimeDb* rtdb=ana->GetRuntimeDb(); + + <CbmMcbm2019CheckDtInDetDataMember> = (<ClassPointer>*) + (rtdb->getContainer("<ContainerName>")); + */ +} + +// ---- Init ---------------------------------------------------------- +InitStatus CbmMcbm2019CheckDtInDet::Init() +{ + + // Get a handle from the IO manager + FairRootManager* ioman = FairRootManager::Instance(); + + // Digi manager + fDigiMan = CbmDigiManager::Instance(); + fDigiMan->UseMuchBeamTimeDigi(); + fDigiMan->Init(); + + // T0 is not included in DigiManager; have to take care here + // Try to find a vector branch for the digi + fT0DigiVector = ioman->InitObjectAs<std::vector<CbmTofDigi> const *>("T0Digi"); + if ( ! fT0DigiVector ) { + LOG(info) << "No T0 digi vector found; trying TClonesArray"; + if ( std::is_convertible<TObject*, CbmTofDigi*>::value ) { + fT0DigiArray = dynamic_cast<TClonesArray*>(ioman->GetObject("T0Digi")); + if ( ! fT0DigiArray ) LOG(info) << "No T0 digi input found."; + } //? CbmTofDigi derives from TObject + } //? No vector for T0 digis + + if ( ! fDigiMan->IsPresent( ECbmModuleId::kSts) ) { + LOG(info) << "No STS digi input found."; + } + + if ( ! fDigiMan->IsPresent( ECbmModuleId::kMuch) ) { + LOG(info) << "No MUCH digi input found."; + } + + if ( ! fDigiMan->IsPresent( ECbmModuleId::kTrd) ) { + LOG(info) << "No TRD digi input found."; + } // if ( ! fDigiMan->IsPresent( ECbmModuleId::kTrd) ) + else { + /// The TRD digi time is relative to the TS start, so we need the metadata to offset it + fTimeSliceMetaDataArray = dynamic_cast<TClonesArray*>(ioman->GetObject("TimesliceMetaData")); + if ( ! fTimeSliceMetaDataArray ) + LOG(fatal) << "No TS metadata input found while TRD needs it."; + } // else of if ( ! fDigiMan->IsPresent( ECbmModuleId::kTrd) ) + + if ( ! fDigiMan->IsPresent( ECbmModuleId::kTof) ) { + LOG(info) << "No TOF digi input found."; + } + + if ( ! fDigiMan->IsPresent( ECbmModuleId::kRich) ) { + LOG(info) << "No RICH digi input found."; + } + + if ( ! fDigiMan->IsPresent( ECbmModuleId::kPsd) ) { + LOG(info) << "No PSD digi input found."; + } + + CreateHistos(); + + return kSUCCESS; +} + +void CbmMcbm2019CheckDtInDet::CreateHistos() +{ + /// Logarithmic bining for self time comparison + uint32_t iNbBinsLog = 0; + /// Parameters are NbDecadesLog, NbStepsDecade, NbSubStepsInStep + std::vector<double> dBinsLogVector = GenerateLogBinArray( 9, 9, 1, iNbBinsLog ); + double* dBinsLog = dBinsLogVector.data(); +// double * dBinsLog = GenerateLogBinArray( 9, 9, 1, iNbBinsLog ); + + /// Proportion of hits with same time + // T0 vs. T0 + fT0T0SameTime = new TH1F("fT0T0SameTime","Fract. same time T0;Same Time? [];Counts", + 2, -0.5, 1.5); + // sts vs. Sts + fStsStsSameTime = new TH1F("fStsStsSameTime","Fract. same time Sts;Same Time? [];Counts", + 2, -0.5, 1.5); + // Much vs. Much + fMuchMuchSameTime = new TH1F("fMuchMuchSameTime","Fract. same time Much;Same Time? [];Counts", + 2, -0.5, 1.5); + // Trd vs. Trd + fTrdTrdSameTime = new TH1F("fTrdTrdSameTime","Fract. same time Trd;Same Time? [];Counts", + 2, -0.5, 1.5); + // Tof vs. Tof + fTofTofSameTime = new TH1F("fTofTofSameTime","Fract. same time Tof;Same Time? [];Counts", + 2, -0.5, 1.5); + // Rich vs. Rich + fRichRichSameTime = new TH1F("fRichRichSameTime","Fract. same time Rich;Same Time? [];Counts", + 2, -0.5, 1.5); + // Psd vs. Psd + fPsdPsdSameTime = new TH1F("fPsdPsdSameTime","Fract. same time Psd;Same Time? [];Counts", + 2, -0.5, 1.5); + + /// Per detector + // T0 vs. T0 + fT0T0Diff = new TH1F("fT0T0Diff","T0-T0_prev;time diff [ns];Counts", + 10001, -0.5, 10000.5); + // sts vs. Sts + fStsStsDiff = new TH1F("fStsStsDiff","Sts-Sts_prev;time diff [ns];Counts", + 10001, -0.5, 10000.5); + // Much vs. Much + fMuchMuchDiff = new TH1F("fMuchMuchDiff","Much-Much_prev;time diff [ns];Counts", + 10001, -0.5, 10000.5); + // Trd vs. Trd + fTrdTrdDiff = new TH1F("fTrdTrdDiff","Trd-Trd_prev;time diff [ns];Counts", + 10001, -0.5, 10000.5); + // Tof vs. Tof + fTofTofDiff = new TH1F("fTofTofDiff","Tof-Tof_prev;time diff [ns];Counts", + 10001, -0.5, 10000.5); + // Rich vs. Rich + fRichRichDiff = new TH1F("fRichRichDiff","Rich-Rich_prev;time diff [ns];Counts", + 10001, -0.5, 10000.5); + // Psd vs. Psd + fPsdPsdDiff = new TH1F("fPsdPsdDiff","Psd-Psd_prev;time diff [ns];Counts", + 10001, -0.5, 10000.5); + // T0 vs. T0 + fT0T0DiffLog = new TH1F("fT0T0DiffLog","T0-T0_prev;time diff [ns];Counts", + iNbBinsLog, dBinsLog); + // sts vs. Sts + fStsStsDiffLog = new TH1F("fStsStsDiffLog","Sts-Sts_prev;time diff [ns];Counts", + iNbBinsLog, dBinsLog); + // Much vs. Much + fMuchMuchDiffLog = new TH1F("fMuchMuchDiffLog","Much-Much_prev;time diff [ns];Counts", + iNbBinsLog, dBinsLog); + // Trd vs. Trd + fTrdTrdDiffLog = new TH1F("fTrdTrdDiffLog","Trd-Trd_prev;time diff [ns];Counts", + iNbBinsLog, dBinsLog); + // Tof vs. Tof + fTofTofDiffLog = new TH1F("fTofTofDiffLog","Tof-Tof_prev;time diff [ns];Counts", + iNbBinsLog, dBinsLog); + // Rich vs. Rich + fRichRichDiffLog = new TH1F("fRichRichDiffLog","Rich-Rich_prev;time diff [ns];Counts", + iNbBinsLog, dBinsLog); + // Psd vs. Psd + fPsdPsdDiffLog = new TH1F("fPsdPsdDiffLog","Psd-Psd_prev;time diff [ns];Counts", + iNbBinsLog, dBinsLog); + + /// Per channel + // T0 vs. T0 + fT0T0DiffPerChan = new TH2F("fT0T0DiffPerChan","T0-T0_prev Per Channel;time diff [ns]; Channel [];Counts", + iNbBinsLog, dBinsLog, + fuNbChanT0, 0, fuNbChanT0 ); + // sts vs. Sts + fStsStsDiffPerChan = new TH2F("fStsStsDiffPerChan","Sts-Sts_prev Per Channel;time diff [ns]; Channel [];Counts", + iNbBinsLog, dBinsLog, + fuNbChanSts, 0, fuNbChanSts ); + // Much vs. Much + fMuchMuchDiffPerChan = new TH2F("fMuchMuchDiffPerChan","Much-Much_prev Per Channel;time diff [ns]; Channel [];Counts", + iNbBinsLog, dBinsLog, + fuNbChanMuch, 0, fuNbChanMuch ); + // Trd vs. Trd + fTrdTrdDiffPerChan = new TH2F("fTrdTrdDiffPerChan","Trd-Trd_prev Per Channel;time diff [ns]; Channel [];Counts", + iNbBinsLog, dBinsLog, + fuNbChanTrd, 0, fuNbChanTrd ); + // Tof vs. Tof + fTofTofDiffPerChan = new TH2F("fTofTofDiffPerChan","Tof-Tof_prev Per Channel;time diff [ns]; Channel [];Counts", + iNbBinsLog, dBinsLog, + fuNbChanTof, 0, fuNbChanTof ); + // Rich vs. Rich + fRichRichDiffPerChan = new TH2F("fRichRichDiffPerChan","Rich-Rich_prev Per Channel;time diff [ns]; Channel [];Counts", + iNbBinsLog, dBinsLog, + fuNbChanRich, 0, fuNbChanRich ); + // Psd vs. Psd + fPsdPsdDiffPerChan = new TH2F("fPsdPsdDiffPerChan","Psd-Psd_prev Per Channel;time diff [ns]; Channel [];Counts", + iNbBinsLog, dBinsLog, + fuNbChanPsd, 0, fuNbChanPsd ); + + /// Register the histos in the HTTP server + FairRunOnline* run = FairRunOnline::Instance(); + if (run) { + THttpServer* server = run->GetHttpServer(); + if( nullptr != server ) { + + server->Register( "/Dt", fT0T0Diff ); + server->Register( "/Dt", fStsStsDiff ); + server->Register( "/Dt", fMuchMuchDiff ); + server->Register( "/Dt", fTrdTrdDiff ); + server->Register( "/Dt", fTofTofDiff ); + server->Register( "/Dt", fRichRichDiff ); + server->Register( "/Dt", fPsdPsdDiff ); + + server->Register( "/Dt", fT0T0DiffLog ); + server->Register( "/Dt", fStsStsDiffLog ); + server->Register( "/Dt", fMuchMuchDiffLog ); + server->Register( "/Dt", fTrdTrdDiffLog ); + server->Register( "/Dt", fTofTofDiffLog ); + server->Register( "/Dt", fRichRichDiffLog ); + server->Register( "/Dt", fPsdPsdDiffLog ); + + server->Register( "/DtPerChan", fT0T0DiffPerChan ); + server->Register( "/DtPerChan", fStsStsDiffPerChan ); + server->Register( "/DtPerChan", fMuchMuchDiffPerChan ); + server->Register( "/DtPerChan", fTrdTrdDiffPerChan ); + server->Register( "/DtPerChan", fTofTofDiffPerChan ); + server->Register( "/DtPerChan", fRichRichDiffPerChan ); + server->Register( "/DtPerChan", fPsdPsdDiffPerChan ); + + } + } +} +// ---- ReInit ------------------------------------------------------- +InitStatus CbmMcbm2019CheckDtInDet::ReInit() +{ + return kSUCCESS; +} + +// ---- Exec ---------------------------------------------------------- +void CbmMcbm2019CheckDtInDet::Exec(Option_t* /*option*/) +{ + LOG(debug) << "executing TS " << fNrTs; + + if( 0 < fNrTs && 0 == fNrTs % 1000 ) + LOG(info) << Form( "Processing TS %6d", fNrTs ); + + /// Get nb entries per detector + LOG(debug) <<"Begin"; + Int_t nrT0Digis = 0; + if ( fT0DigiVector ) nrT0Digis = fT0DigiVector->size(); + else if ( fT0DigiArray ) nrT0Digis = fT0DigiArray->GetEntriesFast(); + LOG(debug) << "T0Digis: " << nrT0Digis; + + Int_t nrStsDigis = fDigiMan->GetNofDigis( ECbmModuleId::kSts); + Int_t nrMuchDigis = fDigiMan->GetNofDigis( ECbmModuleId::kMuch); + Int_t nrTrdDigis = fDigiMan->GetNofDigis( ECbmModuleId::kTrd); + Int_t nrTofDigis = fDigiMan->GetNofDigis( ECbmModuleId::kTof); + Int_t nrRichDigis = fDigiMan->GetNofDigis( ECbmModuleId::kRich); + Int_t nrPsdDigis = fDigiMan->GetNofDigis( ECbmModuleId::kPsd); + + /// Check dT in T0 + for( Int_t iT0 = 0; iT0 < nrT0Digis; ++iT0 ) + { + + if (iT0%1000 == 0) LOG(debug) << "Executing entry " << iT0; + + const CbmTofDigi* T0Digi = nullptr; + if ( fT0DigiVector) T0Digi = &(fT0DigiVector->at(iT0)); + else if ( fT0DigiArray ) T0Digi = dynamic_cast<CbmTofDigi*>(fT0DigiArray->At(iT0)); + assert(T0Digi); + + Double_t T0Time = T0Digi->GetTime(); + Int_t T0Address = T0Digi->GetAddress(); + + Double_t T0TimeDiff = T0Time - fPrevTimeT0; + + if( 0 < iT0 ) + { + fT0T0Diff->Fill( T0TimeDiff ); + if( 0 < T0TimeDiff) + { + fT0T0SameTime->Fill( 0 ); + fT0T0DiffLog->Fill( T0TimeDiff ); + } // if( 0 < T0TimeDiff) + else fT0T0SameTime->Fill( 1 ); + } // if( 0 < iT0 ) + + fPrevTimeT0 = T0Time; + } // for( Int_t iT0 = 0; iT0 < nrT0Digis; ++iT0 ) + + /// Check dT in the other channels + FillHistosPerDet<CbmStsDigi>( fStsStsSameTime, fStsStsDiff, fStsStsDiffLog, fStsStsDiffPerChan, ECbmModuleId::kSts ); + FillHistosPerDet<CbmMuchBeamTimeDigi>( fMuchMuchSameTime, fMuchMuchDiff, fMuchMuchDiffLog, fMuchMuchDiffPerChan, ECbmModuleId::kMuch ); + FillHistosPerDet<CbmTrdDigi>( fTrdTrdSameTime, fTrdTrdDiff, fTrdTrdDiffLog, fTrdTrdDiffPerChan, ECbmModuleId::kTrd ); + FillHistosPerDet<CbmTofDigi>( fTofTofSameTime, fTofTofDiff, fTofTofDiffLog, fTofTofDiffPerChan, ECbmModuleId::kTof ); + FillHistosPerDet<CbmRichDigi>( fRichRichSameTime, fRichRichDiff, fRichRichDiffLog, fRichRichDiffPerChan, ECbmModuleId::kRich ); + FillHistosPerDet<CbmPsdDigi>( fPsdPsdSameTime, fPsdPsdDiff, fPsdPsdDiffLog, fPsdPsdDiffPerChan, ECbmModuleId::kPsd ); + + fNrTs++; + + if( 0 < fNrTs && 0 == fNrTs % 90000 ) + WriteHistos(); +} + + +template <class Digi> void CbmMcbm2019CheckDtInDet::FillHistosPerDet( + TH1* histoSameTime, TH1* histoDt, + TH1* histoDtLog, TH2* histoDtPerChan, + ECbmModuleId iDetId + ) +{ + UInt_t uNrDigis = fDigiMan->GetNofDigis(iDetId); + + Double_t dPrevTime = -1; + + for( UInt_t uDigiIdx = 0; uDigiIdx < uNrDigis; ++uDigiIdx ) + { + const Digi* digi = fDigiMan->Get<Digi>( uDigiIdx ); + + Double_t dNewDigiTime = digi->GetTime(); + + if( 0 < uDigiIdx ) + { + Double_t dTimeDiff = dNewDigiTime - dPrevTime; + + histoDt->Fill( dTimeDiff ); + if( 0 < dTimeDiff ) + { + histoSameTime->Fill( 0 ); + histoDtLog->Fill( dTimeDiff ); + } // if( 0 < dTimeDiff ) + else histoSameTime->Fill( 1 ); + } // +/* + /// Fill histos per Channel + switch( iDetId ) + { + case kSts: ///< Silicon Tracking System + { + const CbmStsDigi* stsDigi; + try { + stsDigi = + boost::any_cast<const CbmStsDigi*>( digi ); + } catch( ... ) { + LOG( fatal ) << "Failed boost any_cast in CbmMcbm2019CheckPulser::FillSystemOffsetHistos for a digi of type " + << Digi::GetClassName(); + } // try/catch + assert(stsDigi); + UInt_t uAddr = stsDigi->GetAddress(); + UInt_t uChan = stsDigi->GetChannel(); + + break; + } // case kSts: + case kMuch: ///< Muon detection system + { + const CbmMuchBeamTimeDigi* muchDigi; + try { + muchDigi = + boost::any_cast<const CbmMuchBeamTimeDigi*>( digi ); + } catch( ... ) { + LOG( fatal ) << "Failed boost any_cast in CbmMcbm2019CheckPulser::FillSystemOffsetHistos for a digi of type " + << Digi::GetClassName(); + } // try/catch + assert(muchDigi); + UInt_t uAsic = muchDigi->GetNxId(); + UInt_t uChan = muchDigi->GetNxCh(); + + break; + } // case kMuch: + case kTrd: ///< Time-of-flight Detector + { + UInt_t uAddr = digi->GetAddress(); + break; + } // case kTrd: + case kTof: ///< Time-of-flight Detector + { + + break; + } // case kTof: + case kRich: ///< Ring-Imaging Cherenkov Detector + { + + break; + } // case kRich: + case kPsd: ///< Projectile spectator detector + { + UInt_t uAddr = digi->GetAddress(); + + break; + } // case kPsd: + default: + return 0; + } // switch( iDetId ) +*/ + dPrevTime = dNewDigiTime; + } // for( UInt_t uDigiIdx = 0; uDigiIdx < uNrDigis; ++uDigiIdx ) +} + +// ---- Finish -------------------------------------------------------- +void CbmMcbm2019CheckDtInDet::Finish() +{ + WriteHistos(); + +} + +void CbmMcbm2019CheckDtInDet::WriteHistos() +{ + TFile* old = gFile; + TFile* outfile = TFile::Open(fOutFileName,"RECREATE"); + + fT0T0SameTime->Write(); + fStsStsSameTime->Write(); + fMuchMuchSameTime->Write(); + fTrdTrdSameTime->Write(); + fTofTofSameTime->Write(); + fRichRichSameTime->Write(); + fPsdPsdSameTime->Write(); + + fT0T0Diff->Write(); + fStsStsDiff->Write(); + fMuchMuchDiff->Write(); + fTrdTrdDiff->Write(); + fTofTofDiff->Write(); + fRichRichDiff->Write(); + fPsdPsdDiff->Write(); + + fT0T0DiffLog->Write(); + fStsStsDiffLog->Write(); + fMuchMuchDiffLog->Write(); + fTrdTrdDiffLog->Write(); + fTofTofDiffLog->Write(); + fRichRichDiffLog->Write(); + fPsdPsdDiffLog->Write(); + + fT0T0DiffPerChan->Write(); + fStsStsDiffPerChan->Write(); + fMuchMuchDiffPerChan->Write(); + fTrdTrdDiffPerChan->Write(); + fTofTofDiffPerChan->Write(); + fRichRichDiffPerChan->Write(); + fPsdPsdDiffPerChan->Write(); + + outfile->Close(); + delete outfile; + + gFile = old; +} + +ClassImp(CbmMcbm2019CheckDtInDet) diff --git a/fles/mcbm2018/tasks/CbmMcbm2019CheckDtInDet.h b/fles/mcbm2018/tasks/CbmMcbm2019CheckDtInDet.h new file mode 100644 index 0000000000000000000000000000000000000000..b30c98b5534f864eb0fb1df6e638946282c80101 --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019CheckDtInDet.h @@ -0,0 +1,170 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#ifndef CBMMCBM2019CHECKDTINDET_H +#define CBMMCBM2019CHECKDTINDET_H + +/// CBMROOT headers +#include "CbmDefs.h" +#include "TimesliceMetaData.h" +#include "CbmTofDigi.h" + +/// FAIRROOT headers +#include "FairTask.h" + +/// External headers (ROOT, boost, ...) +#include "TString.h" + +/// C/C++ headers +#include <vector> + +class TClonesArray; +class TH1; +class TH2; +class TProfile; +class CbmDigiManager; + +class CbmMcbm2019CheckDtInDet : public FairTask +{ + public: + + CbmMcbm2019CheckDtInDet(); + + CbmMcbm2019CheckDtInDet(const CbmMcbm2019CheckDtInDet&) = delete; + CbmMcbm2019CheckDtInDet operator=(const CbmMcbm2019CheckDtInDet&) = delete; + + /** Constructor with parameters (Optional) **/ + // CbmMcbm2019CheckDtInDet(Int_t verbose); + + + /** Destructor **/ + ~CbmMcbm2019CheckDtInDet(); + + + /** Initiliazation of task at the beginning of a run **/ + virtual InitStatus Init(); + + /** ReInitiliazation of task when the runID changes **/ + virtual InitStatus ReInit(); + + + /** Executed for each event. **/ + virtual void Exec(Option_t*); + + /** Load the parameter container from the runtime database **/ + virtual void SetParContainers(); + + /** Finish task called at the end of the run **/ + virtual void Finish(); + + void SetNbChanT0(Int_t val = 8) + { fuNbChanT0 = val; } + + void SetNbChanSts(Int_t val = 5120) + { fuNbChanSts = val; } + + void SetNbChanMuch(Int_t val = 1000) + { fuNbChanMuch = val; } + + void SetNbChanTrd(Int_t val = 1000) + { fuNbChanTrd = val; } + + void SetNbChanTof(Int_t val = 1000) + { fuNbChanTof = val; } + + void SetNbChanRich(Int_t val = 1000) + { fuNbChanRich = val; } + + void SetNbChanPsd(Int_t val = 1000) + { fuNbChanPsd = val; } + + inline void SetOutFilename( TString sNameIn ) { fOutFileName = sNameIn; } + + private: + template <class Digi> void FillHistosPerDet( TH1* histoSameTime, TH1* histoDt, + TH1* histoDtLog ,TH2* histoDtPerChan, + ECbmModuleId iDetId = ECbmModuleId::kLastModule + ); + void CreateHistos(); + void WriteHistos(); + + + /** Digi data **/ + CbmDigiManager* fDigiMan = nullptr; //! + const std::vector<CbmTofDigi>* fT0DigiVector = nullptr; //! + TClonesArray* fT0DigiArray = nullptr; //! + TClonesArray* fTimeSliceMetaDataArray = nullptr; //! + const TimesliceMetaData* pTsMetaData = nullptr; + + /// Constants + static const UInt_t kuNbChanSMX = 128; + static const UInt_t kuMaxNbStsDpbs = 2; + static const UInt_t kuMaxNbMuchDpbs = 6; + static const UInt_t kuMaxNbMuchAsics = 36; + static const UInt_t kuDefaultAddress = 0xFFFFFFFF; + static const UInt_t kuMaxChannelSts = 3000; + + /// Variables to store the previous digi time + Double_t fPrevTimeT0 = 0.; + Double_t fPrevTimeSts = 0.; + Double_t fPrevTimeMuch = 0.; + Double_t fPrevTimeTrd = 0.; + Double_t fPrevTimeTof = 0.; + Double_t fPrevTimeRich = 0.; + Double_t fPrevTimePsd = 0.; + + /// User settings: Data correction parameters + UInt_t fuNbChanT0 = 8; + UInt_t fuNbChanSts = 5120; + UInt_t fuNbChanMuch = 5120; + UInt_t fuNbChanTrd = 5120; + UInt_t fuNbChanTof = 5120; + UInt_t fuNbChanRich = 5120; + UInt_t fuNbChanPsd = 5120; + // + Int_t fNrTs = 0; + + TH1* fT0T0SameTime = nullptr; + TH1* fStsStsSameTime = nullptr; + TH1* fMuchMuchSameTime = nullptr; + TH1* fTrdTrdSameTime = nullptr; + TH1* fTofTofSameTime = nullptr; + TH1* fRichRichSameTime = nullptr; + TH1* fPsdPsdSameTime = nullptr; + + TH1* fT0T0Diff = nullptr; + TH1* fStsStsDiff = nullptr; + TH1* fMuchMuchDiff = nullptr; + TH1* fTrdTrdDiff = nullptr; + TH1* fTofTofDiff = nullptr; + TH1* fRichRichDiff = nullptr; + TH1* fPsdPsdDiff = nullptr; + + TH1* fT0T0DiffLog = nullptr; + TH1* fStsStsDiffLog = nullptr; + TH1* fMuchMuchDiffLog = nullptr; + TH1* fTrdTrdDiffLog = nullptr; + TH1* fTofTofDiffLog = nullptr; + TH1* fRichRichDiffLog = nullptr; + TH1* fPsdPsdDiffLog = nullptr; + + + TH2* fT0T0DiffPerChan = nullptr; + TH2* fStsStsDiffPerChan = nullptr; + TH2* fMuchMuchDiffPerChan = nullptr; + TH2* fTrdTrdDiffPerChan = nullptr; + TH2* fTofTofDiffPerChan = nullptr; + TH2* fRichRichDiffPerChan = nullptr; + TH2* fPsdPsdDiffPerChan = nullptr; + + + TString fOutFileName{"data/HistosDtInDet.root"}; + + ClassDef(CbmMcbm2019CheckDtInDet,1); +}; + +#endif // CBMMCBM2019CHECKDTINDET_H diff --git a/fles/mcbm2018/tasks/CbmMcbm2019CheckTimingPairs.cxx b/fles/mcbm2018/tasks/CbmMcbm2019CheckTimingPairs.cxx new file mode 100644 index 0000000000000000000000000000000000000000..3bae8409b67795cb0606cc742d7c60bee3bd8a35 --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019CheckTimingPairs.cxx @@ -0,0 +1,647 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#include "CbmMcbm2019CheckTimingPairs.h" + +#include "CbmDigiManager.h" +#include "CbmFlesHistosTools.h" + +#include "FairLogger.h" +#include "FairRootManager.h" +#include "FairRunOnline.h" + +#include "TClonesArray.h" +#include "TH1.h" +#include "TH2.h" +#include "TProfile.h" +#include "THttpServer.h" +#include <TFile.h> + + +#include <iomanip> +#include <iostream> +#include <type_traits> +using std::fixed; +using std::setprecision; + +// ---- Default constructor ------------------------------------------- +CbmMcbm2019CheckTimingPairs::CbmMcbm2019CheckTimingPairs() + : FairTask("CbmMcbm2019CheckTimingPairs") +{ +} + +// ---- Destructor ---------------------------------------------------- +CbmMcbm2019CheckTimingPairs::~CbmMcbm2019CheckTimingPairs() +{ +} + +// ---- Initialisation ---------------------------------------------- +void CbmMcbm2019CheckTimingPairs::SetParContainers() +{ + // Load all necessary parameter containers from the runtime data base + /* + FairRunAna* ana = FairRunAna::Instance(); + FairRuntimeDb* rtdb=ana->GetRuntimeDb(); + + <CbmMcbm2019CheckTimingPairsDataMember> = (<ClassPointer>*) + (rtdb->getContainer("<ContainerName>")); + */ +} + +// ---- Init ---------------------------------------------------------- +InitStatus CbmMcbm2019CheckTimingPairs::Init() +{ + + // Get a handle from the IO manager + FairRootManager* ioman = FairRootManager::Instance(); + + // Digi manager + fDigiMan = CbmDigiManager::Instance(); + fDigiMan->UseMuchBeamTimeDigi(); + fDigiMan->Init(); + + // T0 is not included in DigiManager; have to take care here + // Try to find a vector branch for the digi + fT0DigiVector = ioman->InitObjectAs<std::vector<CbmTofDigi> const *>("T0Digi"); + if ( ! fT0DigiVector ) { + LOG(info) << "No T0 digi vector found; trying TClonesArray"; + if ( std::is_convertible<TObject*, CbmTofDigi*>::value ) { + fT0DigiArray = dynamic_cast<TClonesArray*>(ioman->GetObject("T0Digi")); + if ( ! fT0DigiArray ) LOG(info) << "No T0 digi input found."; + } //? CbmTofDigi derives from TObject + } //? No vector for T0 digis + + if ( ! fDigiMan->IsPresent(ECbmModuleId::kSts) ) { + LOG(info) << "No STS digi input found."; + } + + if ( ! fDigiMan->IsPresent(ECbmModuleId::kMuch) ) { + LOG(info) << "No MUCH digi input found."; + } + + if ( ! fDigiMan->IsPresent(ECbmModuleId::kTrd) ) { + LOG(info) << "No TRD digi input found."; + } // if ( ! fDigiMan->IsPresent(ECbmModuleId::kTrd) ) + else { + } // else of if ( ! fDigiMan->IsPresent(ECbmModuleId::kTrd) ) + + if ( ! fDigiMan->IsPresent(ECbmModuleId::kTof) ) { + LOG(info) << "No TOF digi input found."; + } + + if ( ! fDigiMan->IsPresent(ECbmModuleId::kRich) ) { + LOG(info) << "No RICH digi input found."; + } + + if ( ! fDigiMan->IsPresent(ECbmModuleId::kPsd) ) { + LOG(info) << "No PSD digi input found."; + } + + /// Access the TS metadata to know TS start tim + fTimeSliceMetaDataArray = dynamic_cast<TClonesArray*>(ioman->GetObject("TimesliceMetaData")); + if ( ! fTimeSliceMetaDataArray ) + LOG(fatal) << "No TS metadata input found"; + + CreateHistos(); + + return kSUCCESS; +} + +void CbmMcbm2019CheckTimingPairs::CreateHistos() +{ + + Double_t dHistLim = kdDefaultTimeWin + kdClockCycle / 2.0; + for( UInt_t uDetA = 0; uDetA < fvsDetectors.size(); ++uDetA ) + { + for( UInt_t uDetB = uDetA; uDetB < fvsDetectors.size(); ++uDetB ) + { + std::string sName = Form( "hDt%s_Vs_%s", fvsDetectors[ uDetA ].c_str(), fvsDetectors[ uDetB ].c_str() ); + std::string sTitle = Form( "Time diff to T0 for %s VS for %s; Dt %s [ns]; dT %s [ns]; Possible pairs []", + fvsDetectors[ uDetA ].c_str(), fvsDetectors[ uDetB ].c_str(), + fvsDetectors[ uDetA ].c_str(), fvsDetectors[ uDetB ].c_str() ); + fhDtADtB.push_back( new TH2D( sName.c_str(), sTitle.c_str(), + kuNbBinsDefault + 1, -dHistLim, dHistLim, + kuNbBinsDefault + 1, -dHistLim, dHistLim ) ); + } // for( UInt_t uDetB = uDetA; uDetB < fvsDetectors.size(); ++uDetB ) + } // for( UInt_t uDetA = 0; uDetA < fvsDetectors.size(); ++uDetA ) + + + /// Register the histos in the HTTP server + FairRunOnline* run = FairRunOnline::Instance(); + if( run ) + { + THttpServer* server = run->GetHttpServer(); + if( nullptr != server ) + { + for( UInt_t uPair = 0; uPair < fhDtADtB.size(); ++uPair ) + { + server->Register("/PairTiming", fhDtADtB[ uPair ]); + } // for( UInt_t uPair = 0; uPar < fhDtADtB.size(); ++uPair ) + } // if( nullptr != server ) + } // if( run ) +} +// ---- ReInit ------------------------------------------------------- +InitStatus CbmMcbm2019CheckTimingPairs::ReInit() +{ + return kSUCCESS; +} + +// ---- Exec ---------------------------------------------------------- +void CbmMcbm2019CheckTimingPairs::Exec(Option_t* /*option*/) +{ + LOG(debug) << "executing TS " << fNrTs; + + if( 0 < fNrTs && 0 == fNrTs % 1000 ) + LOG(info) << Form( "Processing TS %6d", fNrTs ); + + /// Fill buffers of hits in correlation to T0 + UInt_t uNbT0Digis = 0; + if ( fT0DigiVector ) uNbT0Digis = fT0DigiVector->size(); + else if ( fT0DigiArray ) uNbT0Digis = fT0DigiArray->GetEntriesFast(); + UInt_t uNbStsDigis = fDigiMan->GetNofDigis(ECbmModuleId::kSts); + UInt_t uNbMuchDigis = fDigiMan->GetNofDigis(ECbmModuleId::kMuch); + UInt_t uNbTrdDigis = fDigiMan->GetNofDigis(ECbmModuleId::kTrd); + UInt_t uNbTofDigis = fDigiMan->GetNofDigis(ECbmModuleId::kTof); + UInt_t uNbRichDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich); + UInt_t uNbPsdDigis = fDigiMan->GetNofDigis(ECbmModuleId::kPsd); + + fuNbDigisWithCoincT0 = 0; + for( UInt_t uT0Digi = 0; uT0Digi < uNbT0Digis; ++uT0Digi ) + { + const CbmTofDigi* pDigiT0 = nullptr; + if ( fT0DigiVector) + pDigiT0 = &( fT0DigiVector->at( uT0Digi ) ); + else if ( fT0DigiArray ) + pDigiT0 = dynamic_cast<CbmTofDigi*>( fT0DigiArray->At( uT0Digi ) ); + assert(pDigiT0); + + UInt_t uChannel = pDigiT0->GetChannel(); + Double_t dTimeT0 = pDigiT0->GetTime(); + + fuNbCoincDigisSts = 0; + fuNbCoincDigisMuch = 0; + fuNbCoincDigisTrd = 0; + fuNbCoincDigisTof = 0; + fuNbCoincDigisRich = 0; + fuNbCoincDigisPsd = 0; + for( UInt_t uIndexDet = 0; uIndexDet < fvsDetectors.size(); ++uIndexDet ) + { + if( "STS" == fvsDetectors[ uIndexDet ] ) + { + fvuPrevT0FirstDigiDet[ uIndexDet ] = FillCorrBuffer< CbmStsDigi >( dTimeT0, + fvuPrevT0FirstDigiDet[ uIndexDet ], + -fdStsTimeWin, fdStsTimeWin, + fvDigisSts, + ECbmModuleId::kSts ); + } // if( "STS" == fvsDetectors[ uIndexDet ] ) + else if( "MUCH" == fvsDetectors[ uIndexDet ] ) + { + fvuPrevT0FirstDigiDet[ uIndexDet ] = FillCorrBuffer< CbmMuchBeamTimeDigi >( dTimeT0, + fvuPrevT0FirstDigiDet[ uIndexDet ], + -fdMuchTimeWin, fdMuchTimeWin, + fvDigisMuch, + ECbmModuleId::kMuch ); + } // else if( "MUCH" == fvsDetectors[ uIndexDet ] ) + else if( "TRD" == fvsDetectors[ uIndexDet ] ) + { + fvuPrevT0FirstDigiDet[ uIndexDet ] = FillCorrBuffer< CbmTrdDigi >( dTimeT0, + fvuPrevT0FirstDigiDet[ uIndexDet ], + -fdTrdTimeWin, fdTrdTimeWin, + fvDigisTrd, + ECbmModuleId::kTrd ); + } // else if( "TRD" == fvsDetectors[ uIndexDet ] ) + else if( "TOF" == fvsDetectors[ uIndexDet ] ) + { + fvuPrevT0FirstDigiDet[ uIndexDet ] = FillCorrBuffer< CbmTofDigi >( dTimeT0, + fvuPrevT0FirstDigiDet[ uIndexDet ], + -fdTofTimeWin, fdTofTimeWin, + fvDigisTof, + ECbmModuleId::kTof ); + } // else if( "TOF" == fvsDetectors[ uIndexDet ] ) + else if( "RICH" == fvsDetectors[ uIndexDet ] ) + { + fvuPrevT0FirstDigiDet[ uIndexDet ] = FillCorrBuffer< CbmRichDigi >( dTimeT0, + fvuPrevT0FirstDigiDet[ uIndexDet ], + -fdRichTimeWin, fdRichTimeWin, + fvDigisRich, + ECbmModuleId::kRich ); + } // else if( "RICH" == fvsDetectors[ uIndexDet ] ) + else if( "PSD" == fvsDetectors[ uIndexDet ] ) + { + fvuPrevT0FirstDigiDet[ uIndexDet ] = FillCorrBuffer< CbmPsdDigi >( dTimeT0, + fvuPrevT0FirstDigiDet[ uIndexDet ], + -fdPsdTimeWin, fdPsdTimeWin, + fvDigisPsd, + ECbmModuleId::kPsd ); + } // else if( "PSD" == fvsDetectors[ uIndexDet ] ) + else LOG( fatal ) << "CbmMcbm2019CheckTimingPairs => Unknown detector"; + } // for( UInt_t uIndexDet = 0; uIndexDet < fvsDetectors.size(); ++uIndexDet ) + + /// Store also the T0 Hit if any STS or MUCH coincidence + if( 0 < fuNbCoincDigisSts || + 0 < fuNbCoincDigisMuch || + 0 < fuNbCoincDigisTrd || + 0 < fuNbCoincDigisTof || + 0 < fuNbCoincDigisRich || + 0 < fuNbCoincDigisPsd ) + { + fvDigisT0.push_back( *pDigiT0 ); + ++fuNbDigisWithCoincT0; + + /// Make sure we keep both vector in sync at the same size + if( 0 == fuNbCoincDigisSts ) + fvDigisSts.push_back( std::vector< CbmStsDigi >() ); + if( 0 == fuNbCoincDigisMuch ) + fvDigisMuch.push_back( std::vector< CbmMuchBeamTimeDigi >() ); + if( 0 == fuNbCoincDigisTrd ) + fvDigisTrd.push_back( std::vector< CbmTrdDigi >() ); + if( 0 == fuNbCoincDigisTof ) + fvDigisTof.push_back( std::vector< CbmTofDigi >() ); + if( 0 == fuNbCoincDigisRich ) + fvDigisRich.push_back( std::vector< CbmRichDigi >() ); + if( 0 == fuNbCoincDigisPsd ) + fvDigisPsd.push_back( std::vector< CbmPsdDigi >() ); + } // if( 0 < uNbCoincDigisSts || 0 < uNbCoincDigisMuch ) + } // for( UInt_t uT0Digi = 0; uT0Digi < uNbT0Digis; ++uT0Digi ) + + /// Fill plots from buffers f correlated hits + for( UInt_t uIndexT0 = 0; uIndexT0 < fvDigisT0.size(); ++uIndexT0 ) + { + UInt_t uHistoIdx = 0; + for( UInt_t uIndexDetA = 0; uIndexDetA < fvsDetectors.size(); ++uIndexDetA ) + { + for( UInt_t uIndexDetB = uIndexDetA; uIndexDetB < fvsDetectors.size(); ++uIndexDetB ) + { + if( "STS" == fvsDetectors[ uIndexDetA ] ) + { + FillHistosInter< CbmStsDigi >( uIndexT0, uIndexDetA, uIndexDetB, + fvDigisSts[ uIndexT0 ], uHistoIdx ); + } // if( "STS" == fvsDetectors[ uIndexDetA ] ) + else if( "MUCH" == fvsDetectors[ uIndexDetA ] ) + { + FillHistosInter< CbmMuchBeamTimeDigi >( uIndexT0, uIndexDetA, uIndexDetB, + fvDigisMuch[ uIndexT0 ], uHistoIdx ); + } // else if( "MUCH" == fvsDetectors[ uIndexDetA ] ) + else if( "TRD" == fvsDetectors[ uIndexDetA ] ) + { + FillHistosInter< CbmTrdDigi >( uIndexT0, uIndexDetA, uIndexDetB, + fvDigisTrd[ uIndexT0 ], uHistoIdx ); + } // else if( "TRD" == fvsDetectors[ uIndexDetA ] ) + else if( "TOF" == fvsDetectors[ uIndexDetA ] ) + { + FillHistosInter< CbmTofDigi >( uIndexT0, uIndexDetA, uIndexDetB, + fvDigisTof[ uIndexT0 ], uHistoIdx ); + } // else if( "TOF" == fvsDetectors[ uIndexDetA ] ) + else if( "RICH" == fvsDetectors[ uIndexDetA ] ) + { + FillHistosInter< CbmRichDigi >( uIndexT0, uIndexDetA, uIndexDetB, + fvDigisRich[ uIndexT0 ], uHistoIdx ); + } // else if( "RICH" == fvsDetectors[ uIndexDetA ] ) + else if( "PSD" == fvsDetectors[ uIndexDetA ] ) + { + FillHistosInter< CbmPsdDigi >( uIndexT0, uIndexDetA, uIndexDetB, + fvDigisPsd[ uIndexT0 ], uHistoIdx ); + } // else if( "PSD" == fvsDetectors[ uIndexDetA ] ) + else LOG( fatal ) << "CbmMcbm2019CheckTimingPairs => Unknown detector"; + + uHistoIdx++; + } // for( UInt_t uIndexDetB = uIndexDetA; uIndexDetB < fvsDetectors.size(); ++uIndexDetB ) + } // for( UInt_t uIndexDetA = 0; uIndexDetA < fvsDetectors.size(); ++uIndexDetA ) + + /// Cleanup buffers + fvDigisSts[ uIndexT0 ].clear(); + fvDigisMuch[ uIndexT0 ].clear(); + fvDigisTrd[ uIndexT0 ].clear(); + fvDigisTof[ uIndexT0 ].clear(); + fvDigisRich[ uIndexT0 ].clear(); + fvDigisPsd[ uIndexT0 ].clear(); + } // for( UInt_t uIndexT0 = 0; uIndexT0 < fvDigisT0.size(); ++uIndexT0 ) + /// Cleanup buffers + fvDigisT0.clear(); + fvDigisSts.clear(); + fvDigisMuch.clear(); + fvDigisTrd.clear(); + fvDigisTof.clear(); + fvDigisRich.clear(); + fvDigisPsd.clear(); + + fNrTs++; + + if( 0 < fNrTs && 0 == fNrTs % 10000 ) + WriteHistos(); +} + +template <class Digi > +UInt_t CbmMcbm2019CheckTimingPairs::FillCorrBuffer( Double_t dTimeT0, UInt_t uIndexStart, + Double_t dWinStartTime, Double_t dWinStopTime, + std::vector< std::vector< Digi > > & vDigi, + ECbmModuleId iDetId ) +{ + + UInt_t nrDigis = fDigiMan->GetNofDigis(iDetId); + UInt_t uFirstDigiInWin = uIndexStart; + + for( UInt_t iDigi = uIndexStart; iDigi < nrDigis; ++iDigi ) + { + const Digi* digi = fDigiMan->Get<Digi>( iDigi ); + + Double_t dTimeDet = digi->GetTime(); + Double_t dTimeDiff = dTimeDet - dTimeT0; + + if( dTimeDiff < dWinStartTime ) + { + uFirstDigiInWin = iDigi; + continue; + } // if( dTimeDiff < dWinStartTime ) + else if( dWinStopTime < dTimeDiff ) + { + break; + } // else if( dWinStopTime < dTimeDiff ) of if( dTimeDiff < dWinStartTime ) + + switch( iDetId ) + { + case ECbmModuleId::kSts: ///< Silicon Tracking System + { + const CbmStsDigi* stsDigi; + try { + stsDigi = boost::any_cast<const CbmStsDigi*>( digi ); + } catch( ... ) { + LOG( fatal ) << "Failed boost any_cast in CbmMcbm2019CheckPulser::FillSystemOffsetHistos for a digi of type " + << Digi::GetClassName(); + } // try/catch + assert(stsDigi); + UInt_t uAddr = stsDigi->GetAddress(); + UInt_t uChan = stsDigi->GetChannel(); + + /// Reject pulser digis + if( ( kuDefaultAddress != fuStsAddress && uAddr == fuStsAddress ) ) + continue; + + /// Concidence candidate, store it! + if( 0 == fuNbCoincDigisSts ) + vDigi.push_back( std::vector< Digi >() ); + + vDigi[ fuNbDigisWithCoincT0 ].push_back( ( *digi ) ); + ++fuNbCoincDigisSts; + + break; + } // case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: ///< Muon detection system + { + const CbmMuchBeamTimeDigi* muchDigi{nullptr}; + try { + muchDigi = + boost::any_cast<const CbmMuchBeamTimeDigi*>( digi ); + } catch( ... ) { + LOG( fatal ) << "Failed boost any_cast in CbmMcbm2019CheckPulser::FillSystemOffsetHistos for a digi of type " + << Digi::GetClassName(); + } // try/catch + assert(muchDigi); + UInt_t uAsic = muchDigi->GetNxId(); + UInt_t uChan = muchDigi->GetNxCh(); + + /// Reject pulser digis + if( ( kuMaxNbMuchAsics != fuMuchAsic && uAsic == fuMuchAsic ) ) + continue; + + /// Concidence candidate, store it! + if( 0 == fuNbCoincDigisMuch ) + vDigi.push_back( std::vector< Digi >() ); + + vDigi[ fuNbDigisWithCoincT0 ].push_back( ( *digi ) ); + ++fuNbCoincDigisMuch; + + break; + } // case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: ///< Time-of-flight Detector + { +/* + UInt_t uAddr = digi->GetAddress(); + + if( ( kuDefaultAddress != fuTrdAddress && uAddr == fuTrdAddress ) ) + continue; + + /// Reject pulser digis + if( fuMinChargePulserTrd < digi->GetCharge() && digi->GetCharge() < fuMaxChargePulserTrd ) + continue; +*/ + /// Concidence candidate, store it! + if( 0 == fuNbCoincDigisTrd ) + vDigi.push_back( std::vector< Digi >() ); + + vDigi[ fuNbDigisWithCoincT0 ].push_back( ( *digi ) ); + ++fuNbCoincDigisTrd; + + break; + } // case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: ///< Time-of-flight Detector + { + /// Reject pulser digis + if( fuMinTotPulserTof < digi->GetCharge() && digi->GetCharge() < fuMaxTotPulserTof ) + continue; + + /// Concidence candidate, store it! + if( 0 == fuNbCoincDigisTof ) + vDigi.push_back( std::vector< Digi >() ); + + vDigi[ fuNbDigisWithCoincT0 ].push_back( ( *digi ) ); + ++fuNbCoincDigisTof; + + break; + } // case ECbmModuleId::kTof: + case ECbmModuleId::kRich: ///< Ring-Imaging Cherenkov Detector + { + /// Reject pulser digis + if( fuMinTotPulserRich < digi->GetCharge() && digi->GetCharge() < fuMaxTotPulserRich ) + continue; + + /// Concidence candidate, store it! + if( 0 == fuNbCoincDigisRich ) + vDigi.push_back( std::vector< Digi >() ); + + vDigi[ fuNbDigisWithCoincT0 ].push_back( ( *digi ) ); + ++fuNbCoincDigisRich; + + break; + } // case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: ///< Projectile spectator detector + { + UInt_t uAddr = digi->GetAddress(); + + /// Reject pulser digis + if( ( kuDefaultAddress != fuPsdAddress && uAddr == fuPsdAddress ) ) + continue; + if( fuMinAdcPulserPsd < digi->GetCharge() && digi->GetCharge() < fuMaxAdcPulserPsd ) + continue; +/* + if( digi->GetAddress() == (9<<10)+8 ) + continue; +*/ + /// Concidence candidate, store it! + if( 0 == fuNbCoincDigisPsd ) + vDigi.push_back( std::vector< Digi >() ); + + vDigi[ fuNbDigisWithCoincT0 ].push_back( ( *digi ) ); + ++fuNbCoincDigisPsd; + + break; + } // case ECbmModuleId::kPsd: + default: + return 0; + } // switch( iDetId ) + } // for( UInt_t iDigi = uIndexStart; iDigi < nrDigis; ++iDigi ) + + return uFirstDigiInWin; +} + +template <class DigiA > +void CbmMcbm2019CheckTimingPairs::FillHistosInter( UInt_t uIndexT0, UInt_t uIndexDetA, UInt_t uIndexDetB, + std::vector< DigiA > & vCorrDigA, + UInt_t uHistoIdx ) +{ + if( "STS" == fvsDetectors[ uIndexDetB ] ) + { + FillHistos< DigiA, CbmStsDigi >( uIndexT0, uIndexDetA, uIndexDetB, + vCorrDigA, fvDigisSts[ uIndexT0 ], + uHistoIdx ); + } // if( "STS" == fvsDetectors[ uIndexDetB ] ) + else if( "MUCH" == fvsDetectors[ uIndexDetB ] ) + { + FillHistos< DigiA, CbmMuchBeamTimeDigi >( uIndexT0, uIndexDetA, uIndexDetB, + vCorrDigA, fvDigisMuch[ uIndexT0 ], + uHistoIdx ); + } // else if( "MUCH" == fvsDetectors[ uIndexDetB ] ) + else if( "TRD" == fvsDetectors[ uIndexDetB ] ) + { + FillHistos< DigiA, CbmTrdDigi >( uIndexT0, uIndexDetA, uIndexDetB, + vCorrDigA, fvDigisTrd[ uIndexT0 ], + uHistoIdx ); + } // else if( "TRD" == fvsDetectors[ uIndexDetB ] ) + else if( "TOF" == fvsDetectors[ uIndexDetB ] ) + { + FillHistos< DigiA, CbmTofDigi >( uIndexT0, uIndexDetA, uIndexDetB, + vCorrDigA, fvDigisTof[ uIndexT0 ], + uHistoIdx ); + } // else if( "TOF" == fvsDetectors[ uIndexDetB ] ) + else if( "RICH" == fvsDetectors[ uIndexDetB ] ) + { + FillHistos< DigiA, CbmRichDigi >( uIndexT0, uIndexDetA, uIndexDetB, + vCorrDigA, fvDigisRich[ uIndexT0 ], + uHistoIdx ); + } // else if( "RICH" == fvsDetectors[ uIndexDetB ] ) + else if( "PSD" == fvsDetectors[ uIndexDetB ] ) + { + FillHistos< DigiA, CbmPsdDigi >( uIndexT0, uIndexDetA, uIndexDetB, + vCorrDigA, fvDigisPsd[ uIndexT0 ], + uHistoIdx ); + } // else if( "PSD" == fvsDetectors[ uIndexDetB ] ) + else LOG( fatal ) << "CbmMcbm2019CheckTimingPairs => Unknown detector"; +} +template <class DigiA, class DigiB> +void CbmMcbm2019CheckTimingPairs::FillHistos( UInt_t uIndexT0, UInt_t uIndexDetA, UInt_t uIndexDetB, + std::vector< DigiA > & vCorrDigA, + std::vector< DigiB > & vCorrDigB, + UInt_t uHistoIdx ) +{ + Double_t dTimeT0 = fvDigisT0[ uIndexT0 ].GetTime(); +/* + std::vector< DigiA > vCorrDigA; + std::vector< DigiB > vCorrDigB; + + /// Get detector A + if( "STS" == fvsDetectors[ uIndexDetA ] ) + { + vCorrDigA = fvDigisSts[ uIndexT0 ]; + } // if( "STS" == fvsDetectors[ uIndexDetA ] ) + else if( "MUCH" == fvsDetectors[ uIndexDetA ] ) + { + vCorrDigA = fvDigisMuch[ uIndexT0 ]; + } // else if( "MUCH" == fvsDetectors[ uIndexDetA ] ) + else if( "TRD" == fvsDetectors[ uIndexDetA ] ) + { + vCorrDigA = fvDigisTrd[ uIndexT0 ]; + } // else if( "TRD" == fvsDetectors[ uIndexDetA ] ) + else if( "TOF" == fvsDetectors[ uIndexDetA ] ) + { + vCorrDigA = fvDigisTof[ uIndexT0 ]; + } // else if( "TOF" == fvsDetectors[ uIndexDetA ] ) + else if( "RICH" == fvsDetectors[ uIndexDetA ] ) + { + vCorrDigA = fvDigisRich[ uIndexT0 ]; + } // else if( "RICH" == fvsDetectors[ uIndexDetA ] ) + else if( "PSD" == fvsDetectors[ uIndexDetA ] ) + { + vCorrDigA = fvDigisPsd[ uIndexT0 ]; + } // else if( "PSD" == fvsDetectors[ uIndexDetA ] ) + else LOG( fatal ) << "CbmMcbm2019CheckTimingPairs => Unknown detector"; + + if( "STS" == fvsDetectors[ uIndexDetB ] ) + { + vCorrDigB = fvDigisSts[ uIndexT0 ]; + } // if( "STS" == fvsDetectors[ uIndexDetB ] ) + else if( "MUCH" == fvsDetectors[ uIndexDetB ] ) + { + vCorrDigB = fvDigisMuch[ uIndexT0 ]; + } // else if( "MUCH" == fvsDetectors[ uIndexDetB ] ) + else if( "TRD" == fvsDetectors[ uIndexDetB ] ) + { + vCorrDigB = fvDigisTrd[ uIndexT0 ]; + } // else if( "TRD" == fvsDetectors[ uIndexDetB ] ) + else if( "TOF" == fvsDetectors[ uIndexDetB ] ) + { + vCorrDigB = fvDigisTof[ uIndexT0 ]; + } // else if( "TOF" == fvsDetectors[ uIndexDetB ] ) + else if( "RICH" == fvsDetectors[ uIndexDetB ] ) + { + vCorrDigB = fvDigisRich[ uIndexT0 ]; + } // else if( "RICH" == fvsDetectors[ uIndexDetB ] ) + else if( "PSD" == fvsDetectors[ uIndexDetB ] ) + { + vCorrDigB = fvDigisPsd[ uIndexT0 ]; + } // else if( "PSD" == fvsDetectors[ uIndexDetB ] ) + else LOG( fatal ) << "CbmMcbm2019CheckTimingPairs => Unknown detector"; +*/ + for( UInt_t uIdxDetA = 0; uIdxDetA < vCorrDigA.size(); ++uIdxDetA ) + { + Double_t dTimeDetA = vCorrDigA[ uIdxDetA ].GetTime(); + Double_t dDtDetA = dTimeDetA - dTimeT0; + + for( UInt_t uIdxDetB = 0; uIdxDetB < vCorrDigB.size(); ++uIdxDetB ) + { + Double_t dTimeDetB = vCorrDigB[ uIdxDetB ].GetTime(); + Double_t dDtDetB = dTimeDetB - dTimeT0; + + fhDtADtB[ uHistoIdx ]->Fill( dDtDetA, dDtDetB ); + } // for( UInt_t uIdxDetB = 0; uIdxDetB < vCoincDigisDetB[ uEvent ].size(); ++vCoincDigisDetB ) + } // for( UInt_t uIdxDetA = 0; uIdxDetA < vCoincDigisDetA[ uEvent ].size(); ++vCoincDigisDetA ) +} + + +// ---- Finish -------------------------------------------------------- +void CbmMcbm2019CheckTimingPairs::Finish() +{ + WriteHistos(); + +} + +void CbmMcbm2019CheckTimingPairs::WriteHistos() +{ + TFile* old = gFile; + TFile* outfile = TFile::Open(fOutFileName,"RECREATE"); + + for( UInt_t uPair = 0; uPair < fhDtADtB.size(); ++uPair ) + { + fhDtADtB[ uPair ]->Write(); + } // for( UInt_t uPair = 0; uPair < fhDtADtB.size(); ++uPair ) + + outfile->Close(); + delete outfile; + + gFile = old; +} + +ClassImp(CbmMcbm2019CheckTimingPairs) diff --git a/fles/mcbm2018/tasks/CbmMcbm2019CheckTimingPairs.h b/fles/mcbm2018/tasks/CbmMcbm2019CheckTimingPairs.h new file mode 100644 index 0000000000000000000000000000000000000000..176043c811aedb248146d1212bc0243118b683a3 --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019CheckTimingPairs.h @@ -0,0 +1,223 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#ifndef CBMMCBM2019CHECKTIMINGPAIRS_H +#define CBMMCBM2019CHECKTIMINGPAIRS_H + +/// CBMROOT headers +#include "CbmDefs.h" +#include "TimesliceMetaData.h" + +#include "CbmStsDigi.h" +#include "CbmMuchBeamTimeDigi.h" +#include "CbmRichDigi.h" +#include "CbmTrdDigi.h" +#include "CbmTofDigi.h" +#include "CbmPsdDigi.h" + +/// FAIRROOT headers +#include "FairTask.h" + +/// External headers (ROOT, boost, ...) +#include "TString.h" + +/// C/C++ headers +#include <vector> + +class TClonesArray; +class TH1; +class TH2; +class TProfile; +class CbmDigiManager; + +class CbmMcbm2019CheckTimingPairs : public FairTask +{ + public: + + CbmMcbm2019CheckTimingPairs(); + + CbmMcbm2019CheckTimingPairs(const CbmMcbm2019CheckTimingPairs&) = delete; + CbmMcbm2019CheckTimingPairs operator=(const CbmMcbm2019CheckTimingPairs&) = delete; + + /** Constructor with parameters (Optional) **/ + // CbmMcbm2019CheckTimingPairs(Int_t verbose); + + + /** Destructor **/ + ~CbmMcbm2019CheckTimingPairs(); + + + /** Initiliazation of task at the beginning of a run **/ + virtual InitStatus Init(); + + /** ReInitiliazation of task when the runID changes **/ + virtual InitStatus ReInit(); + + + /** Executed for each event. **/ + virtual void Exec(Option_t*); + + /** Load the parameter container from the runtime database **/ + virtual void SetParContainers(); + + /** Finish task called at the end of the run **/ + virtual void Finish(); + + void SetStsOffsetSearchRange(Double_t val = 1000) + { fdStsTimeWin = val; } + + void SetMuchOffsetSearchRange(Double_t val = 1000) + { fdMuchTimeWin = val; } + + void SetTrdOffsetSearchRange(Double_t val = 1000) + { fdTrdTimeWin = val; } + + void SetTofOffsetSearchRange(Double_t val = 1000) + { fdTofTimeWin = val; } + + void SetRichOffsetSearchRange(Double_t val = 1000) + { fdRichTimeWin = val; } + + void SetPsdOffsetSearchRange(Double_t val = 1000) + { fdPsdTimeWin = val; } + + inline void SetT0PulserTotLimits( UInt_t uMin, UInt_t uMax ) + { fuMinTotPulserT0 = uMin; fuMaxTotPulserT0 = uMax; } + inline void SetStsPulserAdcLimits( UInt_t uMin, UInt_t uMax ) + { fuMinAdcPulserSts = uMin; fuMaxAdcPulserSts = uMax; } + inline void SetMuchPulserAdcLimits( UInt_t uMin, UInt_t uMax ) + { fuMinAdcPulserMuch = uMin; fuMaxAdcPulserMuch = uMax; } + inline void SetTrdPulserChargeLimits( UInt_t uMin, UInt_t uMax ) + { fuMinChargePulserTrd = uMin; fuMaxChargePulserTrd = uMax; } + inline void SetTofPulserTotLimits( UInt_t uMin, UInt_t uMax ) + { fuMinTotPulserTof = uMin; fuMaxTotPulserTof = uMax; } + inline void SetRichPulserTotLimits( UInt_t uMin, UInt_t uMax ) + { fuMinTotPulserRich = uMin; fuMaxTotPulserRich = uMax; } + inline void SetPsdPulserAdcLimits( UInt_t uMin, UInt_t uMax ) + { fuMinAdcPulserPsd = uMin; fuMaxAdcPulserPsd = uMax; } + + inline void SetOutFilename( TString sNameIn ) { fOutFileName = sNameIn; } + + inline void SetStsAddress( UInt_t uAddress ) { fuStsAddress = uAddress; } + inline void SetMuchAsic( UInt_t uAsic ) { fuMuchAsic = uAsic; } + inline void SetMuchChanRange( UInt_t uFirstChan, UInt_t uLastChan = kuNbChanSMX ) + { fuMuchFirstCha = uFirstChan; fuMuchLastChan = uLastChan; } + inline void SetTrdAddress( UInt_t uAddress ) { fuTrdAddress = uAddress; } + inline void SetPsdAddress( UInt_t uAddress ) { fuPsdAddress = uAddress; } + + private: + + void CreateHistos(); + void WriteHistos(); + + template <class Digi > UInt_t FillCorrBuffer( Double_t dTimeT0, UInt_t uIndexStart, + Double_t dWinStartTime, Double_t dWinStopTime, + std::vector< std::vector< Digi > > & vDigi, + ECbmModuleId iDetId = ECbmModuleId::kLastModule); + template <class DigiA > void FillHistosInter( UInt_t uIndexT0, UInt_t uIndexA, UInt_t uIndexB, + std::vector< DigiA > & vCorrDigA, + UInt_t uHistoIdx ); + template <class DigiA, class DigiB> void FillHistos( UInt_t uIndexT0, UInt_t uIndexA, UInt_t uIndexB, + std::vector< DigiA > & vCorrDigA, + std::vector< DigiB > & vCorrDigB, + UInt_t uHistoIdx ); + + /** Digi data **/ + CbmDigiManager* fDigiMan = nullptr; //! + const std::vector<CbmTofDigi>* fT0DigiVector = nullptr; //! + TClonesArray* fT0DigiArray = nullptr; //! + TClonesArray* fTimeSliceMetaDataArray = nullptr; //! + const TimesliceMetaData* pTsMetaData = nullptr; + + /// Constants + static const UInt_t kuNbChanSMX = 128; + static const UInt_t kuMaxNbStsDpbs = 2; + static const UInt_t kuMaxNbMuchDpbs = 6; + static const UInt_t kuMaxNbMuchAsics = 36; + static const UInt_t kuDefaultAddress = 0xFFFFFFFF; + static const UInt_t kuMaxChannelSts = 3000; + static const UInt_t kuNbBinsDefault = 2000; + static constexpr Double_t kdClockCycle = 3.125; + static constexpr Double_t kdDefaultTimeWin = kdClockCycle * kuNbBinsDefault / 2; + + /// List of detectors + std::vector< std::string > fvsDetectors = { "STS", "MUCH", "TRD", "TOF", "RICH", "PSD" }; //! + UInt_t fuNbDetectors = fvsDetectors.size(); + + /// Variables to store the previous digi time + Double_t fPrevTimeT0 = 0.; + std::vector< Double_t > fvPrevTimeDet = std::vector< Double_t >( fuNbDetectors, 0.0 ); //! + + /// Variables to store the first digi fitting the previous T0 hits + /// => Time-order means the time window for following one can only be in a later digi + std::vector< UInt_t > fvuPrevT0FirstDigiDet = std::vector< UInt_t >( fuNbDetectors, 0 ); //! + + /// Variable to store correlated Digis + std::vector< CbmTofDigi > fvDigisT0 = {}; //! + std::vector< std::vector< CbmStsDigi > > fvDigisSts = {}; //! + std::vector< std::vector< CbmMuchBeamTimeDigi > > fvDigisMuch = {}; //! + std::vector< std::vector< CbmTrdDigi > > fvDigisTrd = {}; //! + std::vector< std::vector< CbmTofDigi > > fvDigisTof = {}; //! + std::vector< std::vector< CbmRichDigi > > fvDigisRich = {}; //! + std::vector< std::vector< CbmPsdDigi > > fvDigisPsd = {}; //! + + /// Variable to store counts of T0 with at least one coincidence + UInt_t fuNbDigisWithCoincT0 = 0; + /// Variable to store counts of T0 with at least one coincidence + UInt_t fuNbCoincDigisSts = 0; + UInt_t fuNbCoincDigisMuch = 0; + UInt_t fuNbCoincDigisTrd = 0; + UInt_t fuNbCoincDigisTof = 0; + UInt_t fuNbCoincDigisRich = 0; + UInt_t fuNbCoincDigisPsd = 0; + + /// User settings: Data correction parameters + /// Charge cut + UInt_t fuMinTotPulserT0 = 182; + UInt_t fuMaxTotPulserT0 = 190; + UInt_t fuMinAdcPulserSts = 90; + UInt_t fuMaxAdcPulserSts = 100; + UInt_t fuMinAdcPulserMuch = 5; + UInt_t fuMaxAdcPulserMuch = 15; + UInt_t fuMinChargePulserTrd = 700000; /// Default: No cut + UInt_t fuMaxChargePulserTrd = 0; /// Default: No cut + UInt_t fuMinTotPulserTof = 182; + UInt_t fuMaxTotPulserTof = 190; + UInt_t fuMinTotPulserRich = 700000; /// Default: No cut + UInt_t fuMaxTotPulserRich = 0; /// Default: No cut + UInt_t fuMinAdcPulserPsd = 700000; /// Default: No cut + UInt_t fuMaxAdcPulserPsd = 0; /// Default: No cut + /// Channel selection + UInt_t fuStsAddress = kuDefaultAddress; + UInt_t fuStsFirstCha = kuMaxChannelSts; + UInt_t fuStsLastChan = kuMaxChannelSts; + UInt_t fuMuchAsic = kuMaxNbMuchAsics; + UInt_t fuMuchFirstCha = kuNbChanSMX; + UInt_t fuMuchLastChan = kuNbChanSMX; + UInt_t fuTrdAddress = kuDefaultAddress; + UInt_t fuPsdAddress = kuDefaultAddress; + + // + Int_t fNrTs = 0; + + Double_t fdStsTimeWin = kdDefaultTimeWin; + Double_t fdMuchTimeWin = kdDefaultTimeWin; + Double_t fdTrdTimeWin = kdDefaultTimeWin; + Double_t fdTofTimeWin = kdDefaultTimeWin; + Double_t fdRichTimeWin = kdDefaultTimeWin; + Double_t fdPsdTimeWin = kdDefaultTimeWin; + + Int_t fBinWidth = 1; + + std::vector< TH1 * > fhDtADtB = {}; //! + + TString fOutFileName{"data/HistosTimingPairs.root"}; + + ClassDef(CbmMcbm2019CheckTimingPairs,1); +}; + +#endif // CBMMCBM2019CHECKTIMINGPAIRS_H diff --git a/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderAlgo.cxx b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderAlgo.cxx new file mode 100644 index 0000000000000000000000000000000000000000..020d3710ea2bd19f1d5e826195bca2b4b317f94f --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderAlgo.cxx @@ -0,0 +1,1206 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#include "CbmMcbm2019TimeWinEventBuilderAlgo.h" + +/// CBM headers +#include "CbmEvent.h" +#include "TimesliceMetaData.h" +#include "CbmStsDigi.h" +#include "CbmMuchBeamTimeDigi.h" +#include "CbmTrdDigi.h" +#include "CbmTofDigi.h" +#include "CbmRichDigi.h" +#include "CbmPsdDigi.h" + +#include "CbmDigiManager.h" + +/// FAIRROOT headers +#include "FairLogger.h" +#include "FairRootManager.h" +#include "FairRunOnline.h" + +/// FAIRSOFT headers (geant, boost, ...) +#include "TClonesArray.h" +#include "TH2.h" +#include "TH1.h" +#include "TCanvas.h" +#include "THttpServer.h" + +/// C/C++ headers + +// ---- Default constructor -------------------------------------------- +CbmMcbm2019TimeWinEventBuilderAlgo::CbmMcbm2019TimeWinEventBuilderAlgo() +{ +} + +// ---- Destructor ----------------------------------------------------- +CbmMcbm2019TimeWinEventBuilderAlgo::~CbmMcbm2019TimeWinEventBuilderAlgo() +{ +} + +// ---- Init ----------------------------------------------------------- +Bool_t CbmMcbm2019TimeWinEventBuilderAlgo::InitAlgo() +{ + LOG(info) << "CbmMcbm2019TimeWinEventBuilderAlgo::InitAlgo => Starting sequence"; + + // Get a handle from the IO manager + FairRootManager* ioman = FairRootManager::Instance(); + + // Get a pointer to the previous already existing data level + fDigiMan = CbmDigiManager::Instance(); + fDigiMan->UseMuchBeamTimeDigi(); + fDigiMan->Init(); + + // T0 is not included in DigiManager + fT0DigiVec = ioman->InitObjectAs<std::vector<CbmTofDigi> const *>("T0Digi"); + if( ! fT0DigiVec ) + { + LOG(info) << "No T0 digi input."; + if( ECbmModuleId::kT0 == fRefDet ) + { + LOG(fatal) << "No digi input for reference detector, stopping there!"; + } // if( ECbmModuleId::kT0 == fRefDet ) + } // if( ! fT0DigiVec ) + + if( ! fDigiMan->IsPresent(ECbmModuleId::kSts) ) + { + LOG(info) << "No STS digi input."; + if( ECbmModuleId::kSts == fRefDet ) + { + LOG(fatal) << "No digi input for reference detector, stopping there!"; + } // if( ECbmModuleId::kT0 == fRefDet ) + } // if( ! fDigiMan->IsPresent(ECbmModuleId::kSts) ) + + if( ! fDigiMan->IsPresent(ECbmModuleId::kMuch) ) + { + LOG(info) << "No MUCH digi input."; + if( ECbmModuleId::kMuch == fRefDet ) + { + LOG(fatal) << "No digi input for reference detector, stopping there!"; + } // if( ECbmModuleId::kT0 == fRefDet ) + } // if( ! fDigiMan->IsPresent(ECbmModuleId::kMuch) ) + + if( ! fDigiMan->IsPresent(ECbmModuleId::kTrd) ) + { + LOG(info) << "No TRD digi input."; + if( ECbmModuleId::kTrd == fRefDet ) + { + LOG(fatal) << "No digi input for reference detector, stopping there!"; + } // if( ECbmModuleId::kT0 == fRefDet ) + } // if( ! fDigiMan->IsPresent(ECbmModuleId::kTrd) ) + + if( ! fDigiMan->IsPresent(ECbmModuleId::kTof) ) + { + LOG(info) << "No TOF digi input."; + if( ECbmModuleId::kTof == fRefDet ) + { + LOG(fatal) << "No digi input for reference detector, stopping there!"; + } // if( ECbmModuleId::kT0 == fRefDet ) + } // if( ! fDigiMan->IsPresent(ECbmModuleId::kTof) ) + + if( ! fDigiMan->IsPresent(ECbmModuleId::kRich) ) + { + LOG(info) << "No RICH digi input."; + if( ECbmModuleId::kRich == fRefDet ) + { + LOG(fatal) << "No digi input for reference detector, stopping there!"; + } // if( ECbmModuleId::kT0 == fRefDet ) + } // if( ! fDigiMan->IsPresent(ECbmModuleId::kRich) ) + + if( ! fDigiMan->IsPresent(ECbmModuleId::kPsd) ) + { + LOG(info) << "No PSD digi input."; + if( ECbmModuleId::kPsd == fRefDet ) + { + LOG(fatal) << "No digi input for reference detector, stopping there!"; + } // if( ECbmModuleId::kT0 == fRefDet ) + } // if( ! fDigiMan->IsPresent(ECbmModuleId::kPsd) ) + + /// Access the TS metadata to know TS start tim + fTimeSliceMetaDataArray = dynamic_cast<TClonesArray*>(ioman->GetObject("TimesliceMetaData")); + if ( ! fTimeSliceMetaDataArray ) + LOG(fatal) << "No TS metadata input found"; + + /// Store the time window for the reference detector + switch( fRefDet ) + { + case ECbmModuleId::kSts: + { + fdRefTimeWinBeg = fdStsTimeWinBeg; + fdRefTimeWinEnd = fdStsTimeWinEnd; + break; + } // case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + { + fdRefTimeWinBeg = fdMuchTimeWinBeg; + fdRefTimeWinEnd = fdMuchTimeWinEnd; + break; + } // case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + { + fdRefTimeWinBeg = fdTrdTimeWinBeg; + fdRefTimeWinEnd = fdTrdTimeWinEnd; + break; + } // case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + { + fdRefTimeWinBeg = fdTofTimeWinBeg; + fdRefTimeWinEnd = fdTofTimeWinEnd; + break; + } // case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + { + fdRefTimeWinBeg = fdRichTimeWinBeg; + fdRefTimeWinEnd = fdRichTimeWinEnd; + break; + } // case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + fdRefTimeWinBeg = fdPsdTimeWinBeg; + fdRefTimeWinEnd = fdPsdTimeWinEnd; + break; + } // case ECbmModuleId::kPsd: + case ECbmModuleId::kT0: + { + fdRefTimeWinBeg = fdT0TimeWinBeg; + fdRefTimeWinEnd = fdT0TimeWinEnd; + break; + } // case ECbmModuleId::kT0: + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::Init => " + << "Trying to use unsupported detectore as reference: " + << fRefDet; + break; + } // default: + } // switch( fRefDet ) + + if( fbFillHistos ) + { + CreateHistograms(); + } // if( fbFillHistos ) + + LOG(info) << "CbmMcbm2019TimeWinEventBuilderAlgo::InitAlgo => Done"; + + return kTRUE; +} + +// ---- ProcessTs ------------------------------------------------------ +void CbmMcbm2019TimeWinEventBuilderAlgo::ProcessTs() +{ + LOG_IF(info, fuNrTs%1000 == 0) <<"Begin of TS " << fuNrTs; + + InitTs(); + + BuildEvents(); + + /// Store last event with trigger if not done by other seed + if( nullptr != fCurrentEvent ) + { + /// TODO: store start time of current event ? + // fCurrentEvent->SetStartTime( fPrevTime ); // Replace Seed time with time of first digi in event? + fCurrentEvent->SetEndTime( fdPrevEvtEndTime ); + fEventVector.push_back( fCurrentEvent ); + fuCurEv++; + + /// Prevent building over TS edge + fCurrentEvent = nullptr; + } // if( nullptr != fCurrentEvent ) + + LOG(debug) << "Found " << fEventVector.size() << " triggered events"; + + if( fbFillHistos ) + { + FillHistos(); + } // if( fbFillHistos ) + + fuNrTs++; +} +void CbmMcbm2019TimeWinEventBuilderAlgo::ClearEventVector() +{ + /// Need to delete the object the pointer points to first + int counter = 0; + for( CbmEvent* event: fEventVector ) + { + LOG(debug) << "Event " << counter << " has " << event->GetNofData() << " digis"; + delete event; + counter++; + } // for( CbmEvent* event: fEventVector) + + fEventVector.clear(); +} +// ---- Finish --------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderAlgo::Finish() +{ + LOG(info) << "Total errors: " << fuErrors; +} + +// --------------------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderAlgo::InitTs() +{ + /// Reset TS based variables (analysis per TS = no building over the border) + fuStartIndexT0 = 0; + fuStartIndexSts = 0; + fuStartIndexMuch = 0; + fuStartIndexTrd = 0; + fuStartIndexTof = 0; + fuStartIndexRich = 0; + fuStartIndexPsd = 0; + fuEndIndexT0 = 0; + fuEndIndexSts = 0; + fuEndIndexMuch = 0; + fuEndIndexTrd = 0; + fuEndIndexTof = 0; + fuEndIndexRich = 0; + fuEndIndexPsd = 0; +} + +void CbmMcbm2019TimeWinEventBuilderAlgo::BuildEvents() +{ + /// Call LoopOnSeed with proper template argument + switch( fRefDet ) + { + case ECbmModuleId::kSts: + { + LoopOnSeeds< CbmStsDigi >(); + break; + } // case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + { + LoopOnSeeds< CbmMuchBeamTimeDigi >(); + break; + } // case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + { + LoopOnSeeds< CbmTrdDigi >(); + break; + } // case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + { + LoopOnSeeds< CbmTofDigi >(); + break; + } // case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + { + LoopOnSeeds< CbmRichDigi >(); + break; + } // case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + LoopOnSeeds< CbmPsdDigi >(); + break; + } // case ECbmModuleId::kPsd: + case ECbmModuleId::kT0: + { + LoopOnSeeds< CbmTofDigi >(); + break; + } // case ECbmModuleId::kT0: + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::BuildEvents => " + << "Trying to search event seeds with unsupported det: " + << fRefDet; + break; + } // default: + } // switch( *det ) +} + +template< class DigiSeed > +void CbmMcbm2019TimeWinEventBuilderAlgo::LoopOnSeeds() +{ + pTsMetaData = dynamic_cast<TimesliceMetaData*>(fTimeSliceMetaDataArray->At(0)); + if( nullptr == pTsMetaData ) + LOG(fatal) << Form( "CbmMcbm2019TimeWinEventBuilderAlgo::LoopOnSeeds => No TS metadata found for TS %6u.", fuNrTs ); + + /// Print warning in first TS if time window borders out of potential overlap + if( ( 0.0 < fdEarliestTimeWinBeg && pTsMetaData->GetOverlapDuration() < fdLatestTimeWinEnd ) || + ( pTsMetaData->GetOverlapDuration() < fdWidestTimeWinRange ) ) + { + LOG(warning) << "CbmMcbm2019TimeWinEventBuilderAlgo::LoopOnSeeds => " + << Form( "Event window not fitting in TS overlap, risk of incomplete events: %f %f %f %llu", + fdEarliestTimeWinBeg, fdLatestTimeWinEnd, fdWidestTimeWinRange, pTsMetaData->GetOverlapDuration() ); + } // if end of event window does not fit in overlap for a seed at edge of TS core + + /// Define an acceptance window for the seeds in order to use the overlap + /// part of the TS to avoid incomplete events + Double_t dSeedWindowBeg = pTsMetaData->GetStartTime() + ( 0.0 < fdEarliestTimeWinBeg ? 0.0 : -fdEarliestTimeWinBeg ); + Double_t dSeedWindowEnd = pTsMetaData->GetOverlapStartTime() + ( 0.0 < fdEarliestTimeWinBeg ? 0.0 : -fdEarliestTimeWinBeg ); + if( fbIgnoreTsOverlap ) + { + dSeedWindowBeg = pTsMetaData->GetStartTime(); + dSeedWindowEnd = pTsMetaData->GetOverlapStartTime(); + } // if( fbIgnoreTsOverlap ) + + switch( fRefDet ) + { + case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + UInt_t uNbRefDigis = ( 0 < fDigiMan->GetNofDigis( fRefDet ) ? fDigiMan->GetNofDigis( fRefDet ) : 0 ); + /// Loop on size of vector + for( UInt_t uDigi = 0; uDigi < uNbRefDigis; ++uDigi ) + { + LOG( debug ) << Form( "Checking seed %6u / %6u", uDigi, uNbRefDigis ); + const DigiSeed * pDigi = fDigiMan->Get< DigiSeed >( uDigi ); + /// Check that _entry is not out of range + if( nullptr != pDigi ) + { + Double_t dTime = pDigi->GetTime(); + + /// Check if seed in acceptance window + if( dTime < dSeedWindowBeg ) + { + continue; + } // if( dTime < dSeedWindowBeg ) + else if( dSeedWindowEnd < dTime ) + { + break; + } // else if( dSeedWindowEnd < dTime ) + + /// Check Seed and build event if needed + CheckSeed( dTime, uDigi ); + } // if( nullptr != pDigi ) + } // for( UInt_t uDigi = 0; uDigi < uNbRefDigis; ++uDigi ) + break; + } // Digi containers controlled by DigiManager + case ECbmModuleId::kT0: + { + if ( fT0DigiVec ) + { + /// Loop on size of vector + UInt_t uNbRefDigis = fT0DigiVec->size(); + /// Loop on size of vector + for( UInt_t uDigi = 0; uDigi < uNbRefDigis; ++uDigi ) + { + LOG( debug ) << Form( "Checking seed %6u / %6u", uDigi, uNbRefDigis ); + + Double_t dTime = fT0DigiVec->at( uDigi ).GetTime(); + + /// Check Seed and build event if needed + CheckSeed( dTime, uDigi ); + } // for( UInt_t uDigi = 0; uDigi < uNbRefDigis; ++uDigi ) + } // if ( fT0DigiVec ) + else LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::LoopOnSeeds => " + << "T0 as reference detector but vector not found!"; + break; + } // case ECbmModuleId::kT0 + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::LoopOnSeeds => Unknow reference detector enum! " << fRefDet; + break; + } // default: + } // switch( fRefDet ) +} + +void CbmMcbm2019TimeWinEventBuilderAlgo::CheckSeed( Double_t dSeedTime, UInt_t uSeedDigiIdx ) +{ + /// If previous event valid and event overlap not allowed, check if we are in overlap + /// and react accordingly + if( nullptr != fCurrentEvent && + EOverlapMode::AllowOverlap != fOverMode && + dSeedTime - fdPrevEvtTime < fdWidestTimeWinRange ) + { + /// Within overlap range + switch( fOverMode ) + { + case EOverlapMode::NoOverlap: + { + /// No overlap allowed => reject + LOG( debug1 ) << "Reject seed due to overlap"; + return; + break; + } // case EOverlapMode::NoOverlap: + case EOverlapMode::MergeOverlap: + { + /// Merge overlap mode => do nothing and go on filling current event + break; + } // case EOverlapMode::MergeOverlap: + case EOverlapMode::AllowOverlap: + { + /// Not in Merge overlap mode => should have been catched before, nothing to do + break; + } // case EOverlapMode::AllowOverlap: + } // switch( fOverMode ) + } // if( prev Event exists and overlap not allowed and overlap present ) + else + { + /// Out of overlap range or in overlap allowed mode + /// => store previous event if not empty and create new one + if( nullptr != fCurrentEvent ) + { + /// TODO: store start time of current event ? +// fCurrentEvent->SetStartTime( fPrevTime ); // Replace Seed time with time of first digi in event? + fCurrentEvent->SetEndTime( fdPrevEvtEndTime ); + fEventVector.push_back( fCurrentEvent ); + fuCurEv++; + } // if( nullptr != fCurrentEvent ) + fCurrentEvent = new CbmEvent( fuCurEv, dSeedTime, 0.); + } // else of if( prev Event exists and overlap not allowed and overlap present ) + + /// If window open for reference detector, search for other reference Digis matching it + /// Otherwise only add the current seed + if( fdRefTimeWinBeg < fdRefTimeWinEnd ) + { + switch( fRefDet ) + { + case ECbmModuleId::kSts: + { + SearchMatches< CbmStsDigi >( dSeedTime, fRefDet, fuStartIndexSts ); + break; + } // case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + { + SearchMatches< CbmMuchBeamTimeDigi >( dSeedTime, fRefDet, fuStartIndexMuch ); + break; + } // case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + { + SearchMatches< CbmTrdDigi >( dSeedTime, fRefDet, fuStartIndexTrd ); + break; + } // case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + { + SearchMatches< CbmTofDigi >( dSeedTime, fRefDet, fuStartIndexTof ); + break; + } // case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + { + SearchMatches< CbmRichDigi >( dSeedTime, fRefDet, fuStartIndexRich ); + break; + } // case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + SearchMatches< CbmPsdDigi >( dSeedTime, fRefDet, fuStartIndexPsd ); + break; + } // case ECbmModuleId::kPsd: + case ECbmModuleId::kT0: + { + SearchMatches< CbmTofDigi >( dSeedTime, fRefDet, fuStartIndexT0 ); + break; + } // case ECbmModuleId::kT0: + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::LoopOnSeeds => " + << "Trying to search matches with unsupported det: " + << fRefDet; + break; + } // default: + } // switch( fRefDet ) + + /// Also add the seed if the window starts after the seed + if( 0 < fdRefTimeWinBeg ) + AddDigiToEvent( fRefDet, uSeedDigiIdx ); + } // if( fdRefTimeWinBeg < fdRefTimeWinEnd ) + else + { + AddDigiToEvent( fRefDet, uSeedDigiIdx ); + } // else of if( fdRefTimeWinBeg < fdRefTimeWinEnd ) + + /// Search for matches for each detector in selection list + for( std::vector< ECbmModuleId >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) + { + switch( *det ) + { + case ECbmModuleId::kSts: + { + SearchMatches< CbmStsDigi >( dSeedTime, *det, fuStartIndexSts ); + break; + } // case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + { + SearchMatches< CbmMuchBeamTimeDigi >( dSeedTime, *det, fuStartIndexMuch ); + break; + } // case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + { + SearchMatches< CbmTrdDigi >( dSeedTime, *det, fuStartIndexTrd ); + break; + } // case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + { + SearchMatches< CbmTofDigi >( dSeedTime, *det, fuStartIndexTof ); + break; + } // case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + { + SearchMatches< CbmRichDigi >( dSeedTime, *det, fuStartIndexRich ); + break; + } // case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + SearchMatches< CbmPsdDigi >( dSeedTime, *det, fuStartIndexPsd ); + break; + } // case ECbmModuleId::kPsd: + case ECbmModuleId::kT0: + { + SearchMatches< CbmTofDigi >( dSeedTime, *det, fuStartIndexT0 ); + break; + } // case ECbmModuleId::kT0: + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::LoopOnSeeds => " + << "Trying to search matches with unsupported det: " + << *det; + break; + } // default: + } // switch( *det ) + } // for( std::vector< ECbmModuleId >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) + + /// Check if event is filling trigger conditions and clear it if not + if( HasTrigger( fCurrentEvent ) ) + { + fdPrevEvtTime = dSeedTime; + + /// In case of NoOverlap or MergeOverlap, we can and should start checking teh next window + /// from end of current window in order to save CPU and avoid duplicating + if( EOverlapMode::NoOverlap == fOverMode || + EOverlapMode::MergeOverlap == fOverMode ) + { + fuStartIndexT0 = fuEndIndexT0; + fuStartIndexSts = fuEndIndexSts; + fuStartIndexMuch = fuEndIndexMuch; + fuStartIndexTrd = fuEndIndexTrd; + fuStartIndexTof = fuEndIndexTof; + fuStartIndexRich = fuEndIndexRich; + fuStartIndexPsd = fuEndIndexPsd; + } // If no overlap or merge overlap + } // if( !HasTrigger( fCurrentEvent ) ) + else + { + LOG( debug1 ) << "Reject seed due to Trigger requirements"; + delete fCurrentEvent; + fCurrentEvent = nullptr; /// delete does NOT set a pointer to nullptr... + } // else of if( !HasTrigger( fCurrentEvent ) ) +} + +template < class DigiCheck > +void CbmMcbm2019TimeWinEventBuilderAlgo::SearchMatches( Double_t dSeedTime, ECbmModuleId detMatch, UInt_t & uStartIndex ) +{ + /// This algo relies on time sorted vectors for the selected detectors + UInt_t uLocalIndexStart = uStartIndex; + UInt_t uLocalIndexEnd = uStartIndex; + + /// FIXME: Use method parameters instead to save 1 switch? + Double_t dTimeWinBeg = 0; + Double_t dTimeWinEnd = 0; + switch( detMatch ) + { + case ECbmModuleId::kSts: + { + dTimeWinBeg = fdStsTimeWinBeg; + dTimeWinEnd = fdStsTimeWinEnd; + break; + } // case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + { + dTimeWinBeg = fdMuchTimeWinBeg; + dTimeWinEnd = fdMuchTimeWinEnd; + break; + } // case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + { + dTimeWinBeg = fdTrdTimeWinBeg; + dTimeWinEnd = fdTrdTimeWinEnd; + break; + } // case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + { + dTimeWinBeg = fdTofTimeWinBeg; + dTimeWinEnd = fdTofTimeWinEnd; + break; + } // case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + { + dTimeWinBeg = fdRichTimeWinBeg; + dTimeWinEnd = fdRichTimeWinEnd; + break; + } // case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + dTimeWinBeg = fdPsdTimeWinBeg; + dTimeWinEnd = fdPsdTimeWinEnd; + break; + } // case ECbmModuleId::kPsd: + case ECbmModuleId::kT0: + { + dTimeWinBeg = fdT0TimeWinBeg; + dTimeWinEnd = fdT0TimeWinEnd; + break; + } // case ECbmModuleId::kT0: + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::SearchMatches => " + << "Trying to search matches with unsupported det: " + << detMatch; + break; + } // default: + } // switch( detMatch ) + + /// Check the Digis until out of window + switch( detMatch ) + { + case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + UInt_t uNbSelDigis = ( 0 < fDigiMan->GetNofDigis( detMatch ) ? fDigiMan->GetNofDigis( detMatch ) : 0 ); + /// Loop on size of vector + for( UInt_t uDigi = uStartIndex; uDigi < uNbSelDigis; ++uDigi ) + { + const DigiCheck * pDigi = fDigiMan->Get< DigiCheck >( uDigi ); + /// Check that _entry is not out of range + if( nullptr != pDigi ) + { + Double_t dTime = pDigi->GetTime(); + Double_t dTimeDiff = dTime - dSeedTime; + + LOG( debug4 ) << detMatch << Form( " => Checking match %6u / %6u, dt %f", uDigi, uNbSelDigis, dTimeDiff ); + + /// Check if within time window, update start/stop indices if needed + if( dTimeDiff < dTimeWinBeg ) + { + ++uLocalIndexStart; + continue; + } // if( dTimeDiff < dTimeWinBeg ) + else if( dTimeWinEnd < dTimeDiff ) + { + /// Store as end the first digi out of window to avoid double counting in case of + /// merged overlap event mode + uLocalIndexEnd = uDigi; + break; + } // else if( dTimeWinEnd < dTimeDiff ) of if( dTimeDiff < dTimeWinBeg ) + + AddDigiToEvent( detMatch, uDigi ); + + if( fdPrevEvtEndTime < dTime ) + fdPrevEvtEndTime = dTime; + } // if( nullptr != pDigi ) + } // for( UInt_t uDigi = 0; uDigi < uNbSelDigis; ++uDigi ) + + /// catch the case where we reach the end of the vector before being out of the time window + if( uLocalIndexEnd < uLocalIndexStart ) + uLocalIndexEnd = uNbSelDigis; + + break; + } // Digi containers controlled by DigiManager + case ECbmModuleId::kT0: + { + if ( fT0DigiVec ) + { + /// Loop on size of vector + UInt_t uNbSelDigis = fT0DigiVec->size(); + /// Loop on size of vector + for( UInt_t uDigi = uStartIndex; uDigi < uNbSelDigis; ++uDigi ) + { + Double_t dTime = fT0DigiVec->at( uDigi ).GetTime(); + + Double_t dTimeDiff = dTime - dSeedTime; + + /// Check if within time window, update start/stop indices if needed + if( dTimeDiff < dTimeWinBeg ) + { + ++uLocalIndexStart; + continue; + } // if( dTimeDiff < dTimeWinBeg ) + else if( dTimeWinEnd < dTimeDiff ) + { + /// Store as end the first digi out of window to avoid double counting in case of + /// merged overlap event mode + uLocalIndexEnd = uDigi; + break; + } // else if( dTimeWinEnd < dTimeDiff ) of if( dTimeDiff < dTimeWinBeg ) + + AddDigiToEvent( detMatch, uDigi ); + + if( fdPrevEvtEndTime < dTime ) + fdPrevEvtEndTime = dTime; + } // for( UInt_t uDigi = 0; uDigi < uNbSelDigis; ++uDigi ) + + /// catch the case where we reach the end of the vector before being out of the time window + if( uLocalIndexEnd < uLocalIndexStart ) + uLocalIndexEnd = uNbSelDigis; + } // if ( fT0DigiVec ) + else LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::SearchMatches => " + << "T0 as selection detector but vector not found!"; + + break; + } // case ECbmModuleId::kT0 + default: + { + return; + break; + } // default: + } // switch( detMatch ) + + /// Update the StartIndex and EndIndex for the next event seed + switch( detMatch ) + { + case ECbmModuleId::kSts: + { + fuStartIndexSts = uLocalIndexStart; + fuEndIndexSts = uLocalIndexEnd; + break; + } // case ECbmModuleId::kSts: + case ECbmModuleId::kMuch: + { + fuStartIndexMuch = uLocalIndexStart; + fuEndIndexMuch = uLocalIndexEnd; + break; + } // case ECbmModuleId::kMuch: + case ECbmModuleId::kTrd: + { + fuStartIndexTrd = uLocalIndexStart; + fuEndIndexTrd = uLocalIndexEnd; + break; + } // case ECbmModuleId::kTrd: + case ECbmModuleId::kTof: + { + fuStartIndexTof = uLocalIndexStart; + fuEndIndexTof = uLocalIndexEnd; + break; + } // case ECbmModuleId::kTof: + case ECbmModuleId::kRich: + { + fuStartIndexRich = uLocalIndexStart; + fuEndIndexRich = uLocalIndexEnd; + break; + } // case ECbmModuleId::kRich: + case ECbmModuleId::kPsd: + { + fuStartIndexPsd = uLocalIndexStart; + fuEndIndexPsd = uLocalIndexEnd; + break; + } // case ECbmModuleId::kPsd: + case ECbmModuleId::kT0: + { + fuStartIndexT0 = uLocalIndexStart; + fuEndIndexT0 = uLocalIndexEnd; + break; + } // case ECbmModuleId::kT0: + default: + { + return; + break; + } // default: + } // switch( detMatch ) +} + +void CbmMcbm2019TimeWinEventBuilderAlgo::AddDigiToEvent( ECbmModuleId _system, Int_t _entry ) +{ + // Fill digi index into event + switch (_system) { + case ECbmModuleId::kMvd: fCurrentEvent->AddData( ECbmDataType::kMvdDigi, _entry ); break; + case ECbmModuleId::kSts: fCurrentEvent->AddData( ECbmDataType::kStsDigi, _entry ); break; + case ECbmModuleId::kRich: fCurrentEvent->AddData( ECbmDataType::kRichDigi, _entry ); break; + case ECbmModuleId::kMuch: fCurrentEvent->AddData( ECbmDataType::kMuchDigi, _entry ); break; + case ECbmModuleId::kTrd: fCurrentEvent->AddData( ECbmDataType::kTrdDigi, _entry ); break; + case ECbmModuleId::kTof: fCurrentEvent->AddData( ECbmDataType::kTofDigi, _entry ); break; + case ECbmModuleId::kPsd: fCurrentEvent->AddData( ECbmDataType::kPsdDigi, _entry ); break; + case ECbmModuleId::kT0: fCurrentEvent->AddData( ECbmDataType::kT0Digi, _entry ); break; + default: + break; + } +} + +Bool_t CbmMcbm2019TimeWinEventBuilderAlgo::HasTrigger(CbmEvent* event) +{ + Bool_t hasTrigger{ kTRUE }; + if ( (fT0DigiVec ) && fuTriggerMinT0Digis > 0) { + Int_t iNbDigis = event->GetNofData( ECbmDataType::kT0Digi ); + hasTrigger = hasTrigger && ( -1 != iNbDigis ) && + ( static_cast< UInt_t >( iNbDigis ) >= fuTriggerMinT0Digis ); + if( !hasTrigger ) + { + LOG( debug2 ) << "Event does not have enough T0 digis: " << iNbDigis + << " vs " << fuTriggerMinT0Digis; + return hasTrigger; + } // if( !hasTrigger ) + } + if (fDigiMan->IsPresent(ECbmModuleId::kSts) && fuTriggerMinStsDigis > 0) { + Int_t iNbDigis = event->GetNofData( ECbmDataType::kStsDigi ); + hasTrigger = hasTrigger && ( -1 != iNbDigis ) && + ( static_cast< UInt_t >( iNbDigis ) >= fuTriggerMinStsDigis ); + if( !hasTrigger ) + { + LOG( debug2 ) << "Event does not have enough STS digis: " << iNbDigis + << " vs " << fuTriggerMinStsDigis; + return hasTrigger; + } // if( !hasTrigger ) + } + if (fDigiMan->IsPresent(ECbmModuleId::kMuch) && fuTriggerMinMuchDigis > 0) { + Int_t iNbDigis = event->GetNofData( ECbmDataType::kMuchDigi ); + hasTrigger = hasTrigger && ( -1 != iNbDigis ) && + ( static_cast< UInt_t >( iNbDigis ) >= fuTriggerMinMuchDigis ); + if( !hasTrigger ) + { + LOG( debug2 ) << "Event does not have enough MUCH digis: " << iNbDigis + << " vs " << fuTriggerMinMuchDigis; + return hasTrigger; + } // if( !hasTrigger ) + } + if (fDigiMan->IsPresent(ECbmModuleId::kTrd) && fuTriggerMinTrdDigis > 0) { + Int_t iNbDigis = event->GetNofData( ECbmDataType::kTrdDigi ); + hasTrigger = hasTrigger && ( -1 != iNbDigis ) && + ( static_cast< UInt_t >( iNbDigis ) >= fuTriggerMinTrdDigis ); + if( !hasTrigger ) + { + LOG( debug2 ) << "Event does not have enough TRD digis: " << iNbDigis + << " vs " << fuTriggerMinTrdDigis; + return hasTrigger; + } // if( !hasTrigger ) + } + if (fDigiMan->IsPresent(ECbmModuleId::kTof) && fuTriggerMinTofDigis > 0) { + Int_t iNbDigis = event->GetNofData( ECbmDataType::kTofDigi ); + hasTrigger = hasTrigger && ( -1 != iNbDigis ) && + ( static_cast< UInt_t >( iNbDigis ) >= fuTriggerMinTofDigis); + if( !hasTrigger ) + { + LOG( debug2 ) << "Event does not have enough TOF digis: " << iNbDigis + << " vs " << fuTriggerMinTofDigis; + return hasTrigger; + } // if( !hasTrigger ) + } + if (fDigiMan->IsPresent(ECbmModuleId::kRich) && fuTriggerMinRichDigis > 0) { + Int_t iNbDigis = event->GetNofData( ECbmDataType::kRichDigi ); + hasTrigger = hasTrigger && ( -1 != iNbDigis ) && + ( static_cast< UInt_t >( iNbDigis ) >= fuTriggerMinRichDigis ); + if( !hasTrigger ) + { + LOG( debug2 ) << "Event does not have enough RICH digis: " << iNbDigis + << " vs " << fuTriggerMinRichDigis; + return hasTrigger; + } // if( !hasTrigger ) + } + if (fDigiMan->IsPresent(ECbmModuleId::kPsd) && fuTriggerMinPsdDigis > 0) { + Int_t iNbDigis = event->GetNofData( ECbmDataType::kPsdDigi ); + hasTrigger = hasTrigger && ( -1 != iNbDigis ) && + ( static_cast< UInt_t >( iNbDigis ) >= fuTriggerMinPsdDigis ); + if( !hasTrigger ) + { + LOG( debug2 ) << "Event does not have enough PSD digis: " << iNbDigis + << " vs " << fuTriggerMinPsdDigis; + return hasTrigger; + } // if( !hasTrigger ) + } + + return hasTrigger; +} + +//---------------------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderAlgo::CreateHistograms() +{ + fhEventTime = new TH1F( "hEventTime", "seed time of the events; Seed time [s]; Events", + 60000, 0, 600 ); + fhEventDt = new TH1F( "fhEventDt", "interval in seed time of consecutive events; Seed time [s]; Events", + 2100, -100.5, 1999.5); + fhEventSize = new TH1F( "hEventSize", + "nb of all digis in the event; Nb Digis []; Events []", + 10000, 0, 10000 ); + fhNbDigiPerEvtTime= new TH2I( "hNbDigiPerEvtTime", + "nb of all digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 1000, 0, 10000 ); + + fhNbDigiPerEvtTimeT0 = new TH2I( "hNbDigiPerEvtTimeT0", + "nb of T0 digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 4000, 0, 4000 ); + fhNbDigiPerEvtTimeSts = new TH2I( "hNbDigiPerEvtTimeSts", + "nb of STS digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 4000, 0, 4000 ); + fhNbDigiPerEvtTimeMuch = new TH2I( "hNbDigiPerEvtTimeMuch", + "nb of MUCH digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 4000, 0, 4000 ); + fhNbDigiPerEvtTimeTrd = new TH2I( "hNbDigiPerEvtTimeTrd", + "nb of TRD digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 4000, 0, 4000 ); + fhNbDigiPerEvtTimeTof = new TH2I( "hNbDigiPerEvtTimeTof", + "nb of TOF digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 4000, 0, 4000 ); + fhNbDigiPerEvtTimeRich = new TH2I( "hNbDigiPerEvtTimeRich", + "nb of RICH digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 4000, 0, 4000 ); + fhNbDigiPerEvtTimePsd = new TH2I( "hNbDigiPerEvtTimePsd", + "nb of PSD digis per event vs seed time of the events; Seed time [s]; Nb Digis []; Events []", + 600, 0, 600, + 4000, 0, 4000 ); + + AddHistoToVector( fhEventTime, "evtbuild" ); + AddHistoToVector( fhEventDt, "evtbuild" ); + AddHistoToVector( fhEventSize, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTime, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTimeT0, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTimeSts, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTimeMuch, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTimeTrd, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTimeTof, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTimeRich, "evtbuild" ); + AddHistoToVector( fhNbDigiPerEvtTimePsd, "evtbuild" ); +} +void CbmMcbm2019TimeWinEventBuilderAlgo::FillHistos() +{ + Double_t dPreEvtTime = -1.0; + for( CbmEvent * evt: fEventVector ) + { + fhEventTime->Fill( evt->GetStartTime() * 1e-9 ); + if( 0.0 <= dPreEvtTime ) + { + fhEventDt->Fill( evt->GetStartTime() - dPreEvtTime ); + } // if( 0.0 <= dPreEvtTime ) + fhEventSize->Fill( evt->GetNofData() ); + fhNbDigiPerEvtTime->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData() ); + + fhNbDigiPerEvtTimeT0 ->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData( ECbmDataType::kT0Digi ) ); + fhNbDigiPerEvtTimeSts ->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData( ECbmDataType::kStsDigi ) ); + fhNbDigiPerEvtTimeMuch->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData( ECbmDataType::kMuchDigi ) ); + fhNbDigiPerEvtTimeTrd ->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData( ECbmDataType::kTrdDigi ) ); + fhNbDigiPerEvtTimeTof ->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData( ECbmDataType::kTofDigi ) ); + fhNbDigiPerEvtTimeRich->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData( ECbmDataType::kRichDigi ) ); + fhNbDigiPerEvtTimePsd ->Fill( evt->GetStartTime() * 1e-9, evt->GetNofData( ECbmDataType::kPsdDigi ) ); + + dPreEvtTime = evt->GetStartTime(); + } // for( CbmEvent * evt: fEventVector ) +} +void CbmMcbm2019TimeWinEventBuilderAlgo::ResetHistograms( Bool_t /*bResetTime*/ ) +{ + fhEventTime->Reset(); + fhEventDt ->Reset(); + fhEventSize->Reset(); + fhNbDigiPerEvtTime->Reset(); + fhNbDigiPerEvtTimeT0 ->Reset(); + fhNbDigiPerEvtTimeSts ->Reset(); + fhNbDigiPerEvtTimeMuch->Reset(); + fhNbDigiPerEvtTimeTrd ->Reset(); + fhNbDigiPerEvtTimeTof ->Reset(); + fhNbDigiPerEvtTimeRich->Reset(); + fhNbDigiPerEvtTimePsd ->Reset(); + +/* + if( kTRUE == bResetTime ) + { + /// Also reset the Start time for the evolution plots! + fdStartTime = -1.0; + } // if( kTRUE == bResetTime ) +*/ +} +//---------------------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderAlgo::AddDetector( ECbmModuleId selDet ) +{ + /// FIXME: This is not true in case TOF is used as reference !!!!! + if( fRefDet == selDet ) + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::AddDetector => Cannot add the reference detector as selection detector!"; + } // if( fRefDet == selDet ) + + for( std::vector< ECbmModuleId >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) + { + if( *det == selDet ) + { + LOG( warning ) << "CbmMcbm2019TimeWinEventBuilderAlgo::AddDetector => Doing nothing, selection detector already in list!" << selDet; + return; + } // if( *det == selDet ) + } // for( std::vector< ECbmModuleId >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) + fvDets.push_back( selDet ); +} +void CbmMcbm2019TimeWinEventBuilderAlgo::RemoveDetector( ECbmModuleId selDet ) +{ + for( std::vector< ECbmModuleId >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) + { + if( *det == selDet ) + { + fvDets.erase( det ); + return; + } // if( *det == selDet ) + } // for( std::vector< ECbmModuleId >::iterator det = fvDets.begin(); det != fvDets.end(); ++det ) + LOG( warning ) << "CbmMcbm2019TimeWinEventBuilderAlgo::RemoveDetector => Doing nothing, selection detector not in list!" << selDet; +} +//---------------------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderAlgo::SetTriggerMinNumber( ECbmModuleId selDet, UInt_t uVal ) +{ + /// Store in corresponding members + std::string sDet = ""; + switch( selDet ) + { + case ECbmModuleId::kSts: + { + fuTriggerMinStsDigis = uVal; + sDet = "STS"; + break; + } + case ECbmModuleId::kRich: + { + fuTriggerMinRichDigis = uVal; + sDet = "RICH"; + break; + } + case ECbmModuleId::kMuch: + { + fuTriggerMinMuchDigis = uVal; + sDet = "MUCH"; + break; + } + case ECbmModuleId::kTrd: + { + fuTriggerMinTrdDigis = uVal; + sDet = "TRD"; + break; + } + case ECbmModuleId::kTof: + { + fuTriggerMinTofDigis = uVal; + sDet = "TOF"; + break; + } + case ECbmModuleId::kPsd: + { + fuTriggerMinPsdDigis = uVal; + sDet = "PSD"; + break; + } + case ECbmModuleId::kT0: + { + fuTriggerMinT0Digis = uVal; + sDet = "T0"; + break; + } + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::SetTriggerMinNumber => " + << "Unsupported or unknow detector enum"; + break; + } + } // switch( det ) + + LOG( debug ) << "Set Trigger min limit for " << sDet << " to " << uVal; +} +void CbmMcbm2019TimeWinEventBuilderAlgo::SetTriggerWindow( ECbmModuleId det, Double_t dWinBeg, Double_t dWinEnd ) +{ + /// Check if valid time window: end strictly after beginning + if( dWinEnd <= dWinBeg ) + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::SetTriggerWindow => Invalid time window: [ " + << dWinBeg << ", " << dWinEnd << " ]"; + + std::string sDet = ""; + /// Store in corresponding members + switch( det ) + { + case ECbmModuleId::kSts: + { + fdStsTimeWinBeg = dWinBeg; + fdStsTimeWinEnd = dWinEnd; + sDet = "STS"; + break; + } + case ECbmModuleId::kRich: + { + fdRichTimeWinBeg = dWinBeg; + fdRichTimeWinEnd = dWinEnd; + sDet = "RICH"; + break; + } + case ECbmModuleId::kMuch: + { + fdMuchTimeWinBeg = dWinBeg; + fdMuchTimeWinEnd = dWinEnd; + sDet = "MUCH"; + break; + } + case ECbmModuleId::kTrd: + { + fdTrdTimeWinBeg = dWinBeg; + fdTrdTimeWinEnd = dWinEnd; + sDet = "TRD"; + break; + } + case ECbmModuleId::kTof: + { + fdTofTimeWinBeg = dWinBeg; + fdTofTimeWinEnd = dWinEnd; + sDet = "TOF"; + break; + } + case ECbmModuleId::kPsd: + { + fdPsdTimeWinBeg = dWinBeg; + fdPsdTimeWinEnd = dWinEnd; + sDet = "PSD"; + break; + } + case ECbmModuleId::kT0: + { + fdT0TimeWinBeg = dWinBeg; + fdT0TimeWinEnd = dWinEnd; + sDet = "T0"; + break; + } + default: + { + LOG( fatal ) << "CbmMcbm2019TimeWinEventBuilderAlgo::SetTriggerWindow => " + << "Unsupported or unknow detector enum"; + break; + } + } // switch( det ) + + LOG( debug ) << "Set Trigger window for " << sDet + << " to [ " << dWinBeg << "; " << dWinEnd << " ]"; + + /// Update the variables storing the earliest and latest time window boundaries + UpdateTimeWinBoundariesExtrema(); + /// Update the variable storing the size if widest time window for overlap detection + UpdateWidestTimeWinRange(); +} +void CbmMcbm2019TimeWinEventBuilderAlgo::UpdateTimeWinBoundariesExtrema() +{ + fdEarliestTimeWinBeg = std::min( fdStsTimeWinBeg, std::min( + fdMuchTimeWinBeg, std::min( + fdTrdTimeWinBeg, std::min( + fdTofTimeWinBeg, std::min( + fdRichTimeWinBeg, std::min( + fdPsdTimeWinBeg, + fdT0TimeWinBeg ) ) ) ) ) ); + fdLatestTimeWinEnd = std::max( fdStsTimeWinEnd, std::max( + fdMuchTimeWinEnd, std::max( + fdTrdTimeWinEnd, std::max( + fdTofTimeWinEnd, std::max( + fdRichTimeWinEnd, std::max( + fdPsdTimeWinEnd, + fdT0TimeWinEnd ) ) ) ) ) ); +} +void CbmMcbm2019TimeWinEventBuilderAlgo::UpdateWidestTimeWinRange() +{ + Double_t fdStsTimeWinRange = fdStsTimeWinEnd - fdStsTimeWinBeg; + Double_t fdMuchTimeWinRange = fdMuchTimeWinEnd - fdMuchTimeWinBeg; + Double_t fdTrdTimeWinRange = fdTrdTimeWinEnd - fdTrdTimeWinBeg; + Double_t fdTofTimeWinRange = fdTofTimeWinEnd - fdTofTimeWinBeg; + Double_t fdRichTimeWinRange = fdRichTimeWinEnd - fdRichTimeWinBeg; + Double_t fdPsdTimeWinRange = fdPsdTimeWinEnd - fdPsdTimeWinBeg; + Double_t fdT0TimeWinRange = fdT0TimeWinEnd - fdT0TimeWinBeg; + + fdWidestTimeWinRange = std::max( fdStsTimeWinRange, std::max( + fdMuchTimeWinRange, std::max( + fdTrdTimeWinRange, std::max( + fdTofTimeWinRange, std::max( + fdRichTimeWinRange, std::max( + fdPsdTimeWinRange, + fdT0TimeWinRange ) ) ) ) ) ); +} +//---------------------------------------------------------------------- + +ClassImp(CbmMcbm2019TimeWinEventBuilderAlgo) diff --git a/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderAlgo.h b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderAlgo.h new file mode 100644 index 0000000000000000000000000000000000000000..b9d5fe3fe8fe8e77acae74b55b9901ddb082464d --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderAlgo.h @@ -0,0 +1,223 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#ifndef CBMMCBM2019TIMEWINEVENTBUILDERALGO_H +#define CBMMCBM2019TIMEWINEVENTBUILDERALGO_H + +/// CBM headers +#include "CbmTofDigi.h" + +/// FAIRROOT headers +#include "FairTask.h" + +/// FAIRSOFT headers (geant, boost, ...) + +/// C/C++ headers +#include <array> +#include <map> +#include <set> +#include <tuple> +#include <vector> + +class TimesliceMetaData; +class CbmEvent; +class CbmDigiManager; +class TClonesArray; +class TH1; +class TH2; +class TNamed; +class TCanvas; + +enum class EOverlapMode { NoOverlap, MergeOverlap, AllowOverlap }; + +class CbmMcbm2019TimeWinEventBuilderAlgo +{ + public: + + /** Default constructor **/ + CbmMcbm2019TimeWinEventBuilderAlgo(); + + CbmMcbm2019TimeWinEventBuilderAlgo(const CbmMcbm2019TimeWinEventBuilderAlgo&) = delete; + CbmMcbm2019TimeWinEventBuilderAlgo operator=(const CbmMcbm2019TimeWinEventBuilderAlgo&) = delete; + + /** Destructor **/ + ~CbmMcbm2019TimeWinEventBuilderAlgo(); + + /** Initiliazation at the beginning of a run **/ + Bool_t InitAlgo(); + + /** Executed for each TS. **/ + void ProcessTs(); + + /** Finish called at the end of the run **/ + void Finish(); + + void SetFillHistos( Bool_t var ) { fbFillHistos = var; } + void ResetHistograms( Bool_t bResetTime = kTRUE ); + + void SetReferenceDetector( ECbmModuleId refDet ) { fRefDet = refDet; } + void AddDetector( ECbmModuleId selDet ); + void RemoveDetector( ECbmModuleId selDet ); + + void SetTriggerMinNumber( ECbmModuleId selDet, UInt_t uVal ); + void SetTriggerMinNumberT0( UInt_t uVal ) { fuTriggerMinT0Digis = uVal;} + void SetTriggerMinNumberSts( UInt_t uVal ) { fuTriggerMinStsDigis = uVal;} + void SetTriggerMinNumberMuch( UInt_t uVal ) { fuTriggerMinMuchDigis = uVal;} + void SetTriggerMinNumberTrd( UInt_t uVal ) { fuTriggerMinTrdDigis = uVal;} + void SetTriggerMinNumberTof( UInt_t uVal ) { fuTriggerMinTofDigis = uVal;} + void SetTriggerMinNumberRich( UInt_t uVal ) { fuTriggerMinRichDigis = uVal;} + void SetTriggerMinNumberPsd( UInt_t uVal ) { fuTriggerMinPsdDigis = uVal;} + + void SetTriggerWindow( ECbmModuleId det, Double_t dWinBeg, Double_t dWinEnd ); + + /// Control flags + void SetEventOverlapMode( EOverlapMode mode ) { fOverMode = mode; } + void SetIgnoreTsOverlap( Bool_t bFlagIn = kTRUE ) { fbIgnoreTsOverlap = bFlagIn; } + + /// For monitor algos + void AddHistoToVector( TNamed * pointer, std::string sFolder = "" ) + { fvpAllHistoPointers.push_back( std::pair< TNamed *, std::string >( pointer, sFolder ) ); } + std::vector< std::pair< TNamed *, std::string > > GetHistoVector() { return fvpAllHistoPointers; } + void AddCanvasToVector( TCanvas * pointer, std::string sFolder = "" ) + { fvpAllCanvasPointers.push_back( std::pair< TCanvas *, std::string >( pointer, sFolder ) ); } + std::vector< std::pair< TCanvas *, std::string > > GetCanvasVector() { return fvpAllCanvasPointers; } + + /// Data output access + std::vector< CbmEvent * > & GetEventVector(){ return fEventVector; } + void ClearEventVector(); + + private: + /// Internal methods + void InitTs(); + void BuildEvents(); + + void CreateHistograms(); + void FillHistos(); + + template< class DigiSeed > void LoopOnSeeds(); + void CheckSeed( Double_t dSeedTime, UInt_t uSeedDigiIdx ); + template< class DigiCheck > void SearchMatches( Double_t dSeedTime, ECbmModuleId detMatch, UInt_t & uStartIndex ); + void AddDigiToEvent( ECbmModuleId det, Int_t uIdx ); + Bool_t HasTrigger(CbmEvent*); + + void UpdateTimeWinBoundariesExtrema(); + void UpdateWidestTimeWinRange(); + + /// Constants + static constexpr Double_t kdDefaultTimeWinBeg = -100.0; + static constexpr Double_t kdDefaultTimeWinEnd = 100.0; + + /// User parameters + /// Control flags + Bool_t fbIgnoreTsOverlap = kFALSE; //! Ignore data in Overlap part of the TS + Bool_t fbFillHistos{kTRUE}; //! Switch ON/OFF filling of histograms + /// Event building mode and detectors selection + EOverlapMode fOverMode{ EOverlapMode::AllowOverlap }; + ECbmModuleId fRefDet{ ECbmModuleId::kT0 }; + std::vector< ECbmModuleId > fvDets{ ECbmModuleId::kSts, ECbmModuleId::kMuch, ECbmModuleId::kTrd, + ECbmModuleId::kTof, ECbmModuleId::kRich, ECbmModuleId::kPsd }; + /// Event building trigger parameters + /** Minimum number of T0 digis needed to generate a trigger, 0 means don't use T0 for trigger generation **/ + UInt_t fuTriggerMinT0Digis{0}; + /** Minimum number of Sts digis needed to generate a trigger, 0 means don't use Sts for trigger generation **/ + UInt_t fuTriggerMinStsDigis{0}; + /** Minimum number of Much digis needed to generate a trigger, 0 means don't use Much for trigger generation **/ + UInt_t fuTriggerMinMuchDigis{0}; + /** Minimum number of Trd digis needed to generate a trigger, 0 means don't use Trd for trigger generation **/ + UInt_t fuTriggerMinTrdDigis{0}; + /** Minimum number of Tof digis needed to generate a trigger, 0 means don't use Tof for trigger generation **/ + UInt_t fuTriggerMinTofDigis{0}; + /** Minimum number of Rich digis needed to generate a trigger, 0 means don't use Rich for trigger generation **/ + UInt_t fuTriggerMinRichDigis{0}; + /** Minimum number of Psd digis needed to generate a trigger, 0 means don't use Psd for trigger generation **/ + UInt_t fuTriggerMinPsdDigis{0}; + /// Trigger Window + Double_t fdRefTimeWinBeg = 0.0; + Double_t fdT0TimeWinBeg = 0.0; + Double_t fdStsTimeWinBeg = kdDefaultTimeWinBeg; + Double_t fdMuchTimeWinBeg = kdDefaultTimeWinBeg; + Double_t fdTrdTimeWinBeg = kdDefaultTimeWinBeg; + Double_t fdTofTimeWinBeg = kdDefaultTimeWinBeg; + Double_t fdRichTimeWinBeg = kdDefaultTimeWinBeg; + Double_t fdPsdTimeWinBeg = kdDefaultTimeWinBeg; + + Double_t fdRefTimeWinEnd = 0.0; + Double_t fdT0TimeWinEnd = 0.0; + Double_t fdStsTimeWinEnd = kdDefaultTimeWinEnd; + Double_t fdMuchTimeWinEnd = kdDefaultTimeWinEnd; + Double_t fdTrdTimeWinEnd = kdDefaultTimeWinEnd; + Double_t fdTofTimeWinEnd = kdDefaultTimeWinEnd; + Double_t fdRichTimeWinEnd = kdDefaultTimeWinEnd; + Double_t fdPsdTimeWinEnd = kdDefaultTimeWinEnd; + + Double_t fdEarliestTimeWinBeg = kdDefaultTimeWinBeg; + Double_t fdLatestTimeWinEnd = kdDefaultTimeWinEnd; + Double_t fdWidestTimeWinRange = kdDefaultTimeWinEnd - kdDefaultTimeWinBeg; + + /// Data input + /// FIXME: usage of CbmDigiManager in FairMq context?!? + /// => Maybe by registering vector (or vector reference) to ioman in Device? + CbmDigiManager* fDigiMan = nullptr; //! + const std::vector< CbmTofDigi > * fT0DigiVec = nullptr; //! + TClonesArray* fTimeSliceMetaDataArray = nullptr; //! + const TimesliceMetaData* pTsMetaData = nullptr; + + /// Data ouptut + CbmEvent* fCurrentEvent{nullptr}; //! pointer to the event which is currently build + std::vector<CbmEvent*> fEventVector; //! vector with all created events + + /// Monitoring histograms + /// => Pointers should be filled with TH1*, TH2*, TProfile*, ... + /// ==> To check if object N is of type T, use "T ObjectPointer = dynamic_cast<T>( fvpAllHistoPointers[N].first );" and check for nullptr + /// ==> To get back the original class name use "fvpAllHistoPointers[N].first->ClassName()" which returns a const char * (e.g. "TH1I") + /// ===> Usage example with feeding a THttpServer: + /// ===> #include "TH2.h" + /// ===> std::string sClassName = vHistos[ uHisto ].first.ClassName(); + /// ===> if( !strncmp( sClassName, "TH1", 3 ) ) + /// ===> server->Register( vHistos[ uHisto ].second.data(), dynamic_cast< TH1 * >(vHistos[ uHisto ].first) ); + /// ===> else if( !strncmp( sClassName, "TH2", 3 ) ) + /// ===> server->Register( vHistos[ uHisto ].second.data(), dynamic_cast< TH2 * >(vHistos[ uHisto ].first) ); + std::vector< std::pair< TNamed *, std::string > > fvpAllHistoPointers; //! Vector of pointers to histograms + optional folder name + std::vector< std::pair< TCanvas *, std::string > > fvpAllCanvasPointers; //! Vector of pointers to canvases + optional folder name + + TH1* fhEventTime{nullptr}; //! histogram with the seed time of the events + TH1* fhEventDt{nullptr}; //! histogram with the interval in seed time of consecutive events + TH1* fhEventSize{nullptr}; //! histogram with the nb of all digis in the event + TH2* fhNbDigiPerEvtTime{nullptr}; //! histogram with the nb of all digis per event vs seed time of the events + TH2* fhNbDigiPerEvtTimeT0{nullptr}; //! histogram with the nb of T0 digis per event vs seed time of the events + TH2* fhNbDigiPerEvtTimeSts{nullptr}; //! histogram with the nb of STS digis per event vs seed time of the events + TH2* fhNbDigiPerEvtTimeMuch{nullptr}; //! histogram with the nb of MUCH digis per event vs seed time of the events + TH2* fhNbDigiPerEvtTimeTrd{nullptr}; //! histogram with the nb of TRD digis per event vs seed time of the events + TH2* fhNbDigiPerEvtTimeTof{nullptr}; //! histogram with the nb of TOF digis per event vs seed time of the events + TH2* fhNbDigiPerEvtTimeRich{nullptr}; //! histogram with the nb of RICH digis per event vs seed time of the events + TH2* fhNbDigiPerEvtTimePsd{nullptr}; //! histogram with the nb of PSD digis per event vs seed time of the events + + /// Internal state variables + UInt_t fuCurEv{0}; //! Event Counter + UInt_t fuErrors{0}; //! Error Counter + UInt_t fuNrTs{0}; //! Timeslice Counter + Double_t fdPrevEvtTime{0.}; //! Save previous time information + Double_t fdPrevEvtEndTime{0.}; //! Save previous event last digi time information + UInt_t fuStartIndexT0{0}; //! Start index for first Digi matching previous seed + UInt_t fuStartIndexSts{0}; //! Start index for first Digi matching previous seed + UInt_t fuStartIndexMuch{0}; //! Start index for first Digi matching previous seed + UInt_t fuStartIndexTrd{0}; //! Start index for first Digi matching previous seed + UInt_t fuStartIndexTof{0}; //! Start index for first Digi matching previous seed + UInt_t fuStartIndexRich{0}; //! Start index for first Digi matching previous seed + UInt_t fuStartIndexPsd{0}; //! Start index for first Digi matching previous seed + UInt_t fuEndIndexT0{0}; //! End index for first Digi not matching previous seed + UInt_t fuEndIndexSts{0}; //! End index for first Digi not matching previous seed + UInt_t fuEndIndexMuch{0}; //! End index for first Digi not matching previous seed + UInt_t fuEndIndexTrd{0}; //! End index for first Digi not matching previous seed + UInt_t fuEndIndexTof{0}; //! End index for first Digi not matching previous seed + UInt_t fuEndIndexRich{0}; //! End index for first Digi not matching previous seed + UInt_t fuEndIndexPsd{0}; //! End index for first Digi not matching previous seed + + ClassDefNV(CbmMcbm2019TimeWinEventBuilderAlgo,1); +}; + +#endif // CBMMCBM2019TIMEWINEVENTBUILDERALGO_H diff --git a/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderTask.cxx b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderTask.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9fb0c767c9fc04da36eba5adad89104fd559614d --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderTask.cxx @@ -0,0 +1,230 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#include "CbmMcbm2019TimeWinEventBuilderTask.h" + +#include "CbmEvent.h" + +#include "FairLogger.h" +#include "FairRootManager.h" +#include "FairRunOnline.h" + +#include "TClonesArray.h" +#include "TH2.h" +#include "TH1.h" +#include "THttpServer.h" +#include <TFile.h> + +// ---- Default constructor ------------------------------------------- +CbmMcbm2019TimeWinEventBuilderTask::CbmMcbm2019TimeWinEventBuilderTask() + : FairTask("CbmMcbm2019TimeWinEventBuilderTask") +{ + /// Create Algo. To be made generic/switchable when more event building algo are available! + fpAlgo = new CbmMcbm2019TimeWinEventBuilderAlgo(); +} + +// ---- Destructor ---------------------------------------------------- +CbmMcbm2019TimeWinEventBuilderTask::~CbmMcbm2019TimeWinEventBuilderTask() +{ + +} + +// ---- Initialisation ---------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderTask::SetParContainers() +{ + /// Nothing to do +} + +// ---- Init ---------------------------------------------------------- +InitStatus CbmMcbm2019TimeWinEventBuilderTask::Init() +{ + /// Get a handle from the IO manager + FairRootManager* ioman = FairRootManager::Instance(); + + /// Register output array (CbmEvent) + fEvents = new TClonesArray("CbmEvent",100); + ioman->Register("CbmEvent", "Cbm Event", fEvents, + IsOutputBranchPersistent("CbmEvent")); + + if ( ! fEvents ) LOG(fatal) << "Output branch was not created"; + + /// Call Algo Init method + if( kTRUE == fpAlgo->InitAlgo() ) + return kSUCCESS; + else return kFATAL; +} + +// ---- ReInit ------------------------------------------------------- +InitStatus CbmMcbm2019TimeWinEventBuilderTask::ReInit() +{ + return kSUCCESS; +} + +// ---- Exec ---------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderTask::Exec(Option_t* /*option*/) +{ + LOG(debug2) << "CbmMcbm2019TimeWinEventBuilderTask::Exec => Starting sequence"; + /// Call Algo ProcessTs method + fpAlgo->ProcessTs(); + + /// Save the resulting vector of events in TClonesArray + FillOutput(); + LOG(debug2) << "CbmMcbm2019TimeWinEventBuilderTask::Exec => Done"; +} + + +// ---- Finish -------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderTask::Finish() +{ + if( fbFillHistos ) + { + SaveHistos(); + } // if( fbFillHistos ) + + /// Call Algo finish method + fpAlgo->Finish(); +} + +//---------------------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderTask::FillOutput() +{ + /// Clear TClonesArray before usage. + fEvents->Delete(); + + /// Get vector reference from algo + std::vector<CbmEvent*> vEvents = fpAlgo->GetEventVector(); + + /// Move CbmEvent from temporary vector to TClonesArray + for( CbmEvent* event: vEvents ) + { + LOG(debug) << "Vector: " << event->ToString(); + new ( (*fEvents)[fEvents->GetEntriesFast()] ) CbmEvent(std::move(*event)); + LOG(debug) << "TClonesArray: " + << static_cast<CbmEvent*>(fEvents->At(fEvents->GetEntriesFast()-1))->ToString(); + } // for( CbmEvent* event: vEvents ) + + /// Clear event vector after usage + fpAlgo->ClearEventVector(); +} +//---------------------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderTask::SaveHistos() +{ + /// Obtain vector of pointers on each histo from the algo (+ optionally desired folder) + std::vector< std::pair< TNamed *, std::string > > vHistos = fpAlgo->GetHistoVector(); + + /// (Re-)Create ROOT file to store the histos + TDirectory * oldDir = NULL; + TFile * histoFile = NULL; + /// Store current directory position to allow restore later + oldDir = gDirectory; + /// open separate histo file in recreate mode + histoFile = new TFile( fsOutFileName , "RECREATE"); + histoFile->cd(); + + /// Save all plots and create folders if needed + for( UInt_t uHisto = 0; uHisto < vHistos.size(); ++uHisto ) + { + /// Make sure we end up in chosen folder + TString sFolder = vHistos[ uHisto ].second.data(); + if( nullptr == gDirectory->Get( sFolder ) ) + gDirectory->mkdir( sFolder ); + gDirectory->cd( sFolder ); + + /// Write plot + vHistos[ uHisto ].first->Write(); + + histoFile->cd(); + } // for( UInt_t uHisto = 0; uHisto < vHistos.size(); ++uHisto ) + + /// Restore original directory position + oldDir->cd(); + histoFile->Close(); +} +//---------------------------------------------------------------------- +void CbmMcbm2019TimeWinEventBuilderTask::SetFillHistos( Bool_t bFlag ) +{ + fbFillHistos = bFlag; + if( nullptr != fpAlgo ) + fpAlgo->SetFillHistos( fbFillHistos ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetOutFilename( TString sNameIn ) +{ + fsOutFileName = sNameIn; +} + +void CbmMcbm2019TimeWinEventBuilderTask::SetReferenceDetector( ECbmModuleId refDet ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetReferenceDetector( refDet ); +} +void CbmMcbm2019TimeWinEventBuilderTask::AddDetector( ECbmModuleId selDet ) +{ + if( nullptr != fpAlgo ) + fpAlgo->AddDetector( selDet ); +} +void CbmMcbm2019TimeWinEventBuilderTask::RemoveDetector( ECbmModuleId selDet ) +{ + if( nullptr != fpAlgo ) + fpAlgo->RemoveDetector( selDet ); +} + +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerMinNumberT0( UInt_t uVal ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerMinNumberT0( uVal ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerMinNumberSts( UInt_t uVal ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerMinNumberSts( uVal ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerMinNumberMuch( UInt_t uVal ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerMinNumberMuch( uVal ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerMinNumberTrd( UInt_t uVal ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerMinNumberTrd( uVal ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerMinNumberTof( UInt_t uVal ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerMinNumberTof( uVal ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerMinNumberRich( UInt_t uVal ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerMinNumberRich( uVal ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerMinNumberPsd( UInt_t uVal ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerMinNumberPsd( uVal ); +} + +void CbmMcbm2019TimeWinEventBuilderTask::SetTriggerWindow( ECbmModuleId det, Double_t dWinBeg, Double_t dWinEnd ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetTriggerWindow( det, dWinBeg, dWinEnd ); +} + +void CbmMcbm2019TimeWinEventBuilderTask::SetEventOverlapMode( EOverlapMode mode ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetEventOverlapMode( mode ); +} +void CbmMcbm2019TimeWinEventBuilderTask::SetIgnoreTsOverlap( Bool_t bFlagIn ) +{ + if( nullptr != fpAlgo ) + fpAlgo->SetIgnoreTsOverlap( bFlagIn ); +} + +//---------------------------------------------------------------------- + +ClassImp(CbmMcbm2019TimeWinEventBuilderTask) diff --git a/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderTask.h b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderTask.h new file mode 100644 index 0000000000000000000000000000000000000000..c2b6b98170ca16dd578349d8eaf0f56ffb245ae3 --- /dev/null +++ b/fles/mcbm2018/tasks/CbmMcbm2019TimeWinEventBuilderTask.h @@ -0,0 +1,99 @@ +/******************************************************************************** + * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH * + * * + * This software is distributed under the terms of the * + * GNU Lesser General Public Licence (LGPL) version 3, * + * copied verbatim in the file "LICENSE" * + ********************************************************************************/ +#ifndef CBMMCBM2019TIMEWINEVENTBUILDERTASK_H +#define CBMMCBM2019TIMEWINEVENTBUILDERTASK_H + +/// CBM headers +#include "CbmMcbm2019TimeWinEventBuilderAlgo.h" + +/// FAIRROOT headers +#include "FairTask.h" + +/// FAIRSOFT headers (geant, boost, ...) + +/// C/C++ headers +#include <array> +#include <map> +#include <set> +#include <tuple> +#include <vector> +#include "CbmTofDigi.h" + +class TClonesArray; + +class CbmMcbm2019TimeWinEventBuilderTask : public FairTask +{ + public: + + /** Default constructor **/ + CbmMcbm2019TimeWinEventBuilderTask(); + + CbmMcbm2019TimeWinEventBuilderTask(const CbmMcbm2019TimeWinEventBuilderTask&) = delete; + CbmMcbm2019TimeWinEventBuilderTask operator=(const CbmMcbm2019TimeWinEventBuilderTask&) = delete; + + /** Constructor with parameters (Optional) **/ + // CbmMcbm2019TimeWinEventBuilderTask(Int_t verbose); + + + /** Destructor **/ + ~CbmMcbm2019TimeWinEventBuilderTask(); + + + /** Initiliazation of task at the beginning of a run **/ + virtual InitStatus Init(); + + /** ReInitiliazation of task when the runID changes **/ + virtual InitStatus ReInit(); + + + /** Executed for each event. **/ + virtual void Exec(Option_t*); + + /** Load the parameter container from the runtime database **/ + virtual void SetParContainers(); + + /** Finish task called at the end of the run **/ + virtual void Finish(); + + void SetFillHistos( Bool_t bFlag = kTRUE ); + void SetOutFilename( TString sNameIn ); + + void SetReferenceDetector( ECbmModuleId refDet ); + void AddDetector( ECbmModuleId selDet ); + void RemoveDetector( ECbmModuleId selDet ); + + void SetTriggerMinNumberT0( UInt_t uVal ); + void SetTriggerMinNumberSts( UInt_t uVal ); + void SetTriggerMinNumberMuch( UInt_t uVal ); + void SetTriggerMinNumberTrd( UInt_t uVal ); + void SetTriggerMinNumberTof( UInt_t uVal ); + void SetTriggerMinNumberRich( UInt_t uVal ); + void SetTriggerMinNumberPsd( UInt_t uVal ); + + void SetTriggerWindow( ECbmModuleId det, Double_t dWinBeg, Double_t dWinEnd ); + + void SetEventOverlapMode( EOverlapMode mode ); + void SetIgnoreTsOverlap( Bool_t bFlagIn ); + + private: + void FillOutput(); + void SaveHistos(); + + CbmMcbm2019TimeWinEventBuilderAlgo * fpAlgo = nullptr; + + TClonesArray* fEvents = nullptr; //! output container of CbmEvents + + Bool_t fbFillHistos{kTRUE}; //! Switch ON/OFF filling of histograms + + /** Name of the histogram output file **/ + TString fsOutFileName{"data/HistosEvtWin.root"}; + + ClassDef(CbmMcbm2019TimeWinEventBuilderTask,1); +}; + +#endif // CBMMCBM2019TIMEWINEVENTBUILDERTASK_H