Skip to content
Snippets Groups Projects
Commit 37012ea6 authored by Pierre-Alain Loizeau's avatar Pierre-Alain Loizeau
Browse files

[mCBM 22] Add macro to make events for single det all hits + chain script + plot macro

parent 95d610a8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
macro/beamtime/mcbm2022/mcbm_event_single_det_hits.C 0 → 100644
+ 202
− 0
View file @ 37012ea6
/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Pierre-Alain Loizeau [committer], Adrian Weber */
// --- Includes needed for IDE
#include <RtypesCore.h>
#include <cstdint>
#include <memory>
#include <string>
#include <vector>
#include <math.h>
#include <stdio.h>
#if !defined(__CLING__)
#include "CbmTrdRawMessageSpadic.h"
#include "CbmTrdSpadic.h"
#include <FairLogger.h>
#include <FairRootFileSink.h>
#include <FairRunOnline.h>
#include <Logger.h>
#include <TStopwatch.h>
#include <TSystem.h>
#endif
/// FIXME: Disable clang formatting to keep easy parameters overview
/* clang-format off */
Bool_t mcbm_event_single_det_hits(std::string infile,
UInt_t uRunId,
std::string sSelDet,
uint32_t uMinDigisHit,
double_t dMaxHitWinNs,
std::int32_t nTimeslices = -1,
bool bDigiEvtsOutput = false,
std::string sOutDir = "data/")
{
/// FIXME: Re-enable clang formatting after parameters initial values setting
/* clang-format on */
// --- Logger settings ----------------------------------------------------
TString logLevel = "INFO";
TString logVerbosity = "LOW";
// ------------------------------------------------------------------------
// ----- Environment --------------------------------------------------
TString myName = "mcbm_event"; // this macro's name for screen output
TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
// ------------------------------------------------------------------------
RawEventBuilderDetector selDet = kRawEventBuilderDetUndef;
if ("T0" == sSelDet) { //
selDet = kRawEventBuilderDetT0;
}
else if ("Sts" == sSelDet) {
selDet = kRawEventBuilderDetSts;
}
else if ("Much" == sSelDet) {
selDet = kRawEventBuilderDetMuch;
}
else if ("Trd1D" == sSelDet) {
selDet = kRawEventBuilderDetTrd;
}
else if ("Trd2D" == sSelDet) {
selDet = kRawEventBuilderDetTrd2D;
}
else if ("Tof" == sSelDet) {
selDet = kRawEventBuilderDetTof;
}
else if ("Rich" == sSelDet) {
selDet = kRawEventBuilderDetRich;
}
else {
std::cout << "-E- " << myName << ": Unknown selected system " << sSelDet << std::endl;
return kFALSE;
}
std::vector<std::string> vInFile = {infile};
// ----- Output filename ----------------------------------------------
std::string suffix = ".events.root";
if (bDigiEvtsOutput) suffix = ".digievents.root";
std::string filename = Form("%d_%s_only", uRunId, selDet.sName.c_str()) + suffix;
std::string outfilename = sOutDir + "/" + filename;
std::cout << "-I- " << myName << ": Output file will be " << outfilename << std::endl;
std::string histosfilename = sOutDir + "/" + filename;
if (bDigiEvtsOutput) { //
histosfilename.replace(histosfilename.find(suffix), suffix.size(), ".digievents.hist.root");
}
else {
histosfilename.replace(histosfilename.find(suffix), suffix.size(), ".events.hist.root");
}
std::cout << "-I- " << myName << ": Histos file will be " << histosfilename << std::endl;
// ------------------------------------------------------------------------
// --------------------event builder---------------------------------------
CbmTaskBuildRawEvents* evBuildRaw = new CbmTaskBuildRawEvents();
evBuildRaw->SetEventOverlapMode(EOverlapModeRaw::NoOverlap);
// Remove all detectors + those not there in 2022
evBuildRaw->RemoveDetector(kRawEventBuilderDetT0);
evBuildRaw->RemoveDetector(kRawEventBuilderDetSts);
evBuildRaw->RemoveDetector(kRawEventBuilderDetMuch);
evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd);
evBuildRaw->RemoveDetector(kRawEventBuilderDetTrd2D);
evBuildRaw->RemoveDetector(kRawEventBuilderDetTof);
evBuildRaw->RemoveDetector(kRawEventBuilderDetRich);
evBuildRaw->RemoveDetector(kRawEventBuilderDetPsd);
// Set single detector as reference detector
evBuildRaw->SetReferenceDetector(selDet);
// void SetTsParameters(double TsStartTime, double TsLength, double TsOverLength):
// => TsStartTime=0, TsLength=256ms in 2021, TsOverLength=TS overlap, not used in mCBM2021
//evBuildRaw->SetTsParameters(0.0, 2.56e8, 0.0);, 0.0);
// void SetTsParameters(double TsStartTime, double TsLength, double TsOverLength):
// => TsStartTime=0, TsLength=128 + 1.28 ms in 2022, TsOverLength=1.28 ms (1MS) in mCBM2022
evBuildRaw->SetTsParameters(0.0, 1.28e8, 1.28e6);
/// FIXME: Disable clang formatting to keep easy parameters overview
/* clang-format off */
evBuildRaw->SetTriggerMinNumber(selDet.detId, uMinDigisHit);
evBuildRaw->SetTriggerMaxNumber(selDet.detId, -1);
evBuildRaw->SetTriggerWindow( selDet.detId, -0.005, dMaxHitWinNs);
/// FIXME: Re-enable clang formatting after parameters initial values setting
/* clang-format on */
// Use standard MUCH digis
evBuildRaw->ChangeMuchBeamtimeDigiFlag();
// Enable DigiEvent output if requested
if (bDigiEvtsOutput) {
evBuildRaw->SetDigiEventOutput();
evBuildRaw->SetDigiEventExclusiveTrdExtraction();
}
evBuildRaw->SetOutFilename(histosfilename);
// evBuildRaw->SetOutputBranchPersistent("CbmEvent", kFALSE);
evBuildRaw->SetWriteHistosToFairSink(kFALSE);
// ------------------------------------------------------------------------
// In general, the following parts need not be touched
// ========================================================================
// ----- Timer --------------------------------------------------------
TStopwatch timer;
timer.Start();
// ------------------------------------------------------------------------
// ----- FairRunAna ---------------------------------------------------
auto run = new FairRunAna();
auto inputSource = new FairFileSource(infile);
run->SetSource(inputSource);
auto sink = new FairRootFileSink(outfilename.data());
run->SetSink(sink);
run->SetRunId(uRunId);
run->AddTask(evBuildRaw);
// ------------------------------------------------------------------------
// ----- Logger settings ----------------------------------------------
FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
// ------------------------------------------------------------------------
// ----- Run initialisation -------------------------------------------
std::cout << std::endl;
std::cout << "-I- " << myName << ": Initialise run" << std::endl;
run->Init();
// ------------------------------------------------------------------------
// ----- Start run ----------------------------------------------------
std::cout << std::endl << std::endl;
if (nTimeslices < 0) {
std::cout << "-I- " << myName << ": Starting run over all timeslices in input" << std::endl;
run->Run(0, -1);
}
else {
std::cout << "-I- " << myName << ": Starting run over " << nTimeslices
<< " timeslices (or less if not enough in input)" << std::endl;
run->Run(0, nTimeslices);
}
// ------------------------------------------------------------------------
// ----- Finish -------------------------------------------------------
timer.Stop();
std::cout << "Macro finished successfully." << std::endl;
std::cout << "After CpuTime = " << timer.CpuTime() << " s RealTime = " << timer.RealTime() << " s." << std::endl;
// ------------------------------------------------------------------------
return kTRUE;
} // End of main macro function
macro/beamtime/mcbm2022/mcbm_hits_det_by_det.sh 0 → 100644
+ 61
− 0
View file @ 37012ea6
#!/bin/bash
NB_TS=1
DIGI_FILE_IN="${VMCWORKDIR}/build/macro/run/data/2570_node8_1_0000.digi.root"
RUN_ID=2570
OUT_DIR="${VMCWORKDIR}/build/macro/run/data/"
if [ $# -gt 1 ]; then
NB_TS=$1
if [ $# -gt 2 ]; then
DIGI_FILE_IN=$2
if [ $# -gt 3 ]; then
RUN_ID=$3
if [ $# -eq 4 ]; then
OUT_DIR=$4
else
echo 'Too many parameters. Only following pattern allowed:'
echo 'mcbm_hits_det_by_det.sh [<Nb TS to process> [<Digi file path> [<Run Id> [<Output dir] ] ] ]'
return -1
fi
fi
fi
fi
# Build events from single detector, one by one
EVENT_MACRO=${VMCWORKDIR}/macro/beamtime/mcbm2022/mcbm_event_single_det_hits.C
root -l -b -q "${EVENT_MACRO}(\"${DIGI_FILE_IN}\", ${RUN_ID}, \"T0\", 1, 5.5, ${NB_TS}, false, \"${OUT_DIR}\")"
root -l -b -q "${EVENT_MACRO}(\"${DIGI_FILE_IN}\", ${RUN_ID}, \"Sts\", 2, 15.5, ${NB_TS}, false, \"${OUT_DIR}\")"
root -l -b -q "${EVENT_MACRO}(\"${DIGI_FILE_IN}\", ${RUN_ID}, \"Much\", 1, 30.5, ${NB_TS}, false, \"${OUT_DIR}\")"
root -l -b -q "${EVENT_MACRO}(\"${DIGI_FILE_IN}\", ${RUN_ID}, \"Trd2D\", 1, 160.5, ${NB_TS}, false, \"${OUT_DIR}\")"
root -l -b -q "${EVENT_MACRO}(\"${DIGI_FILE_IN}\", ${RUN_ID}, \"Trd1D\", 1, 250.5, ${NB_TS}, false, \"${OUT_DIR}\")"
root -l -b -q "${EVENT_MACRO}(\"${DIGI_FILE_IN}\", ${RUN_ID}, \"Tof\", 2, 30.5, ${NB_TS}, false, \"${OUT_DIR}\")"
root -l -b -q "${EVENT_MACRO}(\"${DIGI_FILE_IN}\", ${RUN_ID}, \"Rich\", 1, 30.5, ${NB_TS}, false, \"${OUT_DIR}\")"
# Reco the single detector events files
RECO_MACRO=${VMCWORKDIR}/macro/beamtime/mcbm2022/mcbm_reco.C
# STS TRD TRD2D RICH MUCH T0/TOF
# kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE,
# STS TRD TRD2D RICH MUCH T0/TOF
root -l -b -q "${RECO_MACRO}(${RUN_ID}, ${NB_TS}, \"${OUT_DIR}\", \"${OUT_DIR}\", -1, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, \"${OUT_DIR}/2570_T0_only.events.root\", false, \"${DIGI_FILE_IN}\")"
mv ${OUT_DIR}\reco_event_mcbm_test${RUN_ID}.root ${OUT_DIR}\reco_${RUN_ID}_T0_only.root
root -l -b -q "${RECO_MACRO}(${RUN_ID}, ${NB_TS}, \"${OUT_DIR}\", \"${OUT_DIR}\", -1, kFALSE, kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, \"${OUT_DIR}/2570_Sts_only.events.root\", false, \"${DIGI_FILE_IN}\")"
mv ${OUT_DIR}\reco_event_mcbm_test${RUN_ID}.root ${OUT_DIR}\reco_${RUN_ID}_Sts_only.root
root -l -b -q "${RECO_MACRO}(${RUN_ID}, ${NB_TS}, \"${OUT_DIR}\", \"${OUT_DIR}\", -1, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, \"${OUT_DIR}/2570_Much_only.events.root\", false, \"${DIGI_FILE_IN}\")"
mv ${OUT_DIR}\reco_event_mcbm_test${RUN_ID}.root ${OUT_DIR}\reco_${RUN_ID}_Much_only.root
root -l -b -q "${RECO_MACRO}(${RUN_ID}, ${NB_TS}, \"${OUT_DIR}\", \"${OUT_DIR}\", -1, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, \"${OUT_DIR}/2570_Trd2D_only.events.root\", false, \"${DIGI_FILE_IN}\")"
mv ${OUT_DIR}\reco_event_mcbm_test${RUN_ID}.root ${OUT_DIR}\reco_${RUN_ID}_Trd2D_only.root
root -l -b -q "${RECO_MACRO}(${RUN_ID}, ${NB_TS}, \"${OUT_DIR}\", \"${OUT_DIR}\", -1, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, \"${OUT_DIR}/2570_Trd1D_only.events.root\", false, \"${DIGI_FILE_IN}\")"
mv ${OUT_DIR}\reco_event_mcbm_test${RUN_ID}.root ${OUT_DIR}\reco_${RUN_ID}_Trd1D_only.root
root -l -b -q "${RECO_MACRO}(${RUN_ID}, ${NB_TS}, \"${OUT_DIR}\", \"${OUT_DIR}\", -1, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, \"${OUT_DIR}/2570_Tof_only.events.root\", false, \"${DIGI_FILE_IN}\")"
mv ${OUT_DIR}\reco_event_mcbm_test${RUN_ID}.root ${OUT_DIR}\reco_${RUN_ID}_Tof_only.root
root -l -b -q "${RECO_MACRO}(${RUN_ID}, ${NB_TS}, \"${OUT_DIR}\", \"${OUT_DIR}\", -1, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, \"${OUT_DIR}/2570_Rich_only.events.root\", false, \"${DIGI_FILE_IN}\")"
mv ${OUT_DIR}\reco_event_mcbm_test${RUN_ID}.root ${OUT_DIR}\reco_${RUN_ID}_Rich_only.root
# Merge and plot
PLOT_MACRO=${VMCWORKDIR}/macro/beamtime/mcbm2022/mcbm_plot_det_by_det_hits.C
echo root -l "${PLOT_MACRO}(${RUN_ID}, \"${OUT_DIR}/reco_${RUN_ID}_T0_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Sts_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Much_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Trd2D_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Trd1D_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Tof_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Rich_only.root\")"
root -l "${PLOT_MACRO}(${RUN_ID}, \"${OUT_DIR}/reco_${RUN_ID}_T0_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Sts_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Much_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Trd2D_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Trd1D_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Tof_only.root\", \"${OUT_DIR}/reco_${RUN_ID}_Rich_only.root\")"
macro/beamtime/mcbm2022/mcbm_plot_det_by_det_hits.C 0 → 100644
+ 141
− 0
View file @ 37012ea6
/* Copyright (C) 2023 Facility for Antiproton and Ion Research in Europe, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Pierre-Alain Loizeau [committer] */
void mcbm_plot_det_by_det_hits(uint32_t uRunIndex, std::string sFileBmon, std::string sFileSts, std::string sFileMuch,
std::string sFileTrd2d, std::string sFileTrd1d, std::string sFileTof,
std::string sFileRich)
{
gROOT->cd();
std::vector<CbmTofHit> vBmonHits;
std::vector<CbmStsHit> vStsHits;
std::vector<CbmMuchPixelHit> vMuchHits;
std::vector<CbmTrdHit> vTrd2dHits;
std::vector<CbmTrdHit> vTrd1dHits;
std::vector<CbmTofHit> vTofHits;
std::vector<CbmRichHit> vRichHits;
/// Open all files
TFile* fBmon = new TFile(sFileBmon.c_str(), "READ");
TFile* fSts = new TFile(sFileSts.c_str(), "READ");
TFile* fMuch = new TFile(sFileMuch.c_str(), "READ");
TFile* fTrd2d = new TFile(sFileTrd2d.c_str(), "READ");
TFile* fTrd1d = new TFile(sFileTrd1d.c_str(), "READ");
TFile* fTof = new TFile(sFileTof.c_str(), "READ");
TFile* fRich = new TFile(sFileRich.c_str(), "READ");
gROOT->cd();
/// Trees
TTree* pTreeBmon = dynamic_cast<TTree*>(fBmon->Get("cbmsim"));
TTree* pTreeSts = dynamic_cast<TTree*>(fSts->Get("cbmsim"));
TTree* pTreeMuch = dynamic_cast<TTree*>(fMuch->Get("cbmsim"));
TTree* pTreeTrd2d = dynamic_cast<TTree*>(fTrd1d->Get("cbmsim"));
TTree* pTreeTrd1d = dynamic_cast<TTree*>(fTrd1d->Get("cbmsim"));
TTree* pTreeTof = dynamic_cast<TTree*>(fTof->Get("cbmsim"));
TTree* pTreeRich = dynamic_cast<TTree*>(fRich->Get("cbmsim"));
/// Check that they all contains the same number of TS
bool bTsNbAllOk = true;
uint32_t uNbTs = pTreeBmon->GetEntries();
if (uNbTs != pTreeSts->GetEntries()) {
std::cout << "Different nb TS Bmon and STS: " << uNbTs << " VS " << pTreeBmon->GetEntries() << std::endl;
bTsNbAllOk = false;
}
if (uNbTs != pTreeMuch->GetEntries()) {
std::cout << "Different nb TS Bmon and Much: " << uNbTs << " VS " << pTreeMuch->GetEntries() << std::endl;
bTsNbAllOk = false;
}
if (uNbTs != pTreeTrd2d->GetEntries()) {
std::cout << "Different nb TS Bmon and Trd2d: " << uNbTs << " VS " << pTreeTrd2d->GetEntries() << std::endl;
bTsNbAllOk = false;
}
if (uNbTs != pTreeTrd1d->GetEntries()) {
std::cout << "Different nb TS Bmon and Trd1d: " << uNbTs << " VS " << pTreeTrd1d->GetEntries() << std::endl;
bTsNbAllOk = false;
}
if (uNbTs != pTreeTof->GetEntries()) {
std::cout << "Different nb TS Bmon and Tof: " << uNbTs << " VS " << pTreeTof->GetEntries() << std::endl;
bTsNbAllOk = false;
}
if (uNbTs != pTreeRich->GetEntries()) {
std::cout << "Different nb TS Bmon and Rich: " << uNbTs << " VS " << pTreeRich->GetEntries() << std::endl;
bTsNbAllOk = false;
}
if (!bTsNbAllOk) {
std::cout << "At least one detector with different nb TS, stopping there!" << std::endl;
return;
}
/// Branches
TClonesArray* pArrBmon = new TClonesArray("CbmTofHit");
TClonesArray* pArrSts = new TClonesArray("CbmStsHit");
TClonesArray* pArrMuch = new TClonesArray("CbmMuchPixelHit");
TClonesArray* pArrTrd2D = new TClonesArray("CbmTrdHit");
TClonesArray* pArrTrd1d = new TClonesArray("CbmTrdHit");
TClonesArray* pArrTof = new TClonesArray("CbmTofHit");
TClonesArray* pArrRich = new TClonesArray("CbmRichHit");
pTreeBmon->SetBranchAddress("TofHit", &pArrBmon);
pTreeSts->SetBranchAddress("StsHit", &pArrSts);
// pTreeMuch->SetBranchAddress("MuchHit", &pArrMuch);
pTreeTrd1d->SetBranchAddress("TrdHit", &pArrTrd1d);
pTreeTrd2d->SetBranchAddress("TrdHit", &pArrTrd2D);
pTreeTof->SetBranchAddress("TofHit", &pArrTof);
pTreeRich->SetBranchAddress("RichHit", &pArrRich);
/// Plots
TH1* hSts = new TH1I("Sts", "Sts/Bmon; TS; Sts/Bmon", 100, 0., 100.);
TH1* hMuch = new TH1I("Much", "Much/Bmon; TS; Much/Bmon", 100, 0., 100.);
TH1* hTrd2d = new TH1I("Trd2d", "Trd2d/Bmon; TS; Trd2d/Bmon", 100, 0., 100.);
TH1* hTrd1d = new TH1I("Trd1d", "Trd1d/Bmon; TS; Trd1d/Bmon", 100, 0., 100.);
TH1* hTof = new TH1I("Tof", "Tof/Bmon; TS; Tof/Bmon", 100, 0., 100.);
TH1* hRich = new TH1I("Rich", "Rich/Bmon; TS; Rich/Bmon", 100, 0., 100.);
/// Loop On TS
for (uint32_t uTs = 0; uTs < uNbTs; ++uTs) {
/// Load all entries
pTreeBmon->GetEntry(uTs);
pTreeSts->GetEntry(uTs);
pTreeMuch->GetEntry(uTs);
pTreeTrd1d->GetEntry(uTs);
pTreeTrd2d->GetEntry(uTs);
pTreeTof->GetEntry(uTs);
pTreeRich->GetEntry(uTs);
/// Extract all hits for each detector
/// Loop on hits and plot as function of time
hSts->Fill(uTs, pArrSts->GetEntriesFast() / pArrBmon->GetEntriesFast());
hMuch->Fill(uTs, pArrMuch->GetEntriesFast() / pArrBmon->GetEntriesFast());
hTrd1d->Fill(uTs, pArrTrd1d->GetEntriesFast() / pArrBmon->GetEntriesFast());
hTrd2d->Fill(uTs, pArrTrd2D->GetEntriesFast() / pArrBmon->GetEntriesFast());
hTof->Fill(uTs, pArrTof->GetEntriesFast() / pArrBmon->GetEntriesFast());
hRich->Fill(uTs, pArrRich->GetEntriesFast() / pArrBmon->GetEntriesFast());
}
TCanvas* canv = new TCanvas("canv", "canv");
canv->cd();
THStack* stack = new THStack("stack", "nb rel. Bmon");
hSts->SetLineColor(kBlack);
stack->Add(hSts);
hMuch->SetLineColor(kRed);
stack->Add(hMuch);
hTrd2d->SetLineColor(kMagenta);
stack->Add(hTrd2d);
hTrd1d->SetLineColor(kGreen);
stack->Add(hTrd1d);
hTof->SetLineColor(kBlue);
stack->Add(hTof);
hRich->SetLineColor(kGray);
stack->Add(hRich);
stack->Draw("nostack,hist");
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment