Commit f06a9602 authored by Rafał Lalik's avatar Rafał Lalik
Browse files

STS drift time trend tool

parent 5dce9bbc
add_executable(
sts_calibrator
add_executable(sts_calibrator sts_calibrator.cpp)
target_link_libraries(
sts_calibrator
Forward
${HYDRA2_LIBRARIES}
${ROOT_LIBRARIES}
)
sts_calibrator.cpp
add_executable(sts_drift_time_trend sts_drift_time_trend.cpp)
target_link_libraries(
sts_drift_time_trend
Forward
${HYDRA2_LIBRARIES}
${ROOT_LIBRARIES}
)
target_link_libraries(sts_calibrator Forward ${HYDRA2_LIBRARIES} ${ROOT_LIBRARIES})
install(
TARGETS sts_drift_time_trend
RUNTIME
DESTINATION ${CMAKE_INSTALL_BINDIR}
)
#include "forward_tools.h"
#include "ft_config.h"
#include <hades.h>
#include <hloop.h>
#include <htaskset.h>
#include <hcategory.h>
#include <hcategorymanager.h>
#include <hrecevent.h>
#include <hreconstructor.h>
#include <hruntimedb.h>
//--------category definitions---------
#include <hgeantdef.h>
#include <hparticledef.h>
#include <hpiontrackerdef.h>
#include <hstartdef.h>
//-------------------------------------
//-------objects-----------------------
#include <heventheader.h>
#include <hgeantkine.h>
#include <hparticlecand.h>
#include <hparticlecandsim.h>
#include <hparticleevtinfo.h>
#include <hparticletracksorter.h>
#include <hpiontrackertrack.h>
#include <hstart2hit.h>
#include <hstart2cal.h>
//-------------------------------------
#include <hparticletool.h>
#include <hphysicsconstants.h>
#include <henergylosscorrpar.h>
#include <fwdetdef.h>
#include <hdst.h>
#include <hforwardcand.h>
#include <hfrpccal.h>
#include <hfrpccalpar.h>
#include <hfrpccluster.h>
#include <hfrpcdigipar.h>
#include <hfrpcgeompar.h>
#include <hfrpchit.h>
#include <hfrpcraw.h>
#include <hgeantfrpc.h>
#include <hgeantsts.h>
#include <hgeomcompositevolume.h>
#include <hgeomvector.h>
#include <hgeomvolume.h>
#include <hparasciifileio.h>
#include <hparticlecand.h>
#include <hrpccal.h>
#include <hrpccalpar.h>
#include <hrpccluster.h>
#include <hrpcdigipar.h>
#include <hrpcgeompar.h>
#include <hrpchit.h>
#include <hrpcraw.h>
#include <hspectrometer.h>
#include <hstart2hit.h>
#include <hstscal.h>
#include <hstsgeompar.h>
#include <hstsraw.h>
#include <htofhit.h>
#include <htofraw.h>
#include <rpcdef.h>
#include <tofdef.h>
#include <TCanvas.h>
#include <TColor.h>
#include <TFile.h>
#include <TGaxis.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TH1I.h>
#include <TH2F.h>
#include <TH2I.h>
#include <TH3I.h>
#include <THStack.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TMatrixD.h>
#include <TMultiGraph.h>
#include <TROOT.h>
#include <TStopwatch.h>
#include <TString.h>
#include <TStyle.h>
#include <TTree.h>
#include <TVector3.h>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <getopt.h>
#include <iostream>
#include <regex>
#include <sstream>
#include <string>
#include <vector>
#define PR(x) \
std::cout << "++DEBUG: " << #x << " = |" << x << "| (" << __FILE__ << ", " << __LINE__ << ")\n";
using namespace std;
using namespace ConfigBeamtime;
//##############################################################################################
// Usefull functions and definitions
//##############################################################################################
//----------------------------------------------------------------------------------------------
// Defining parameters
//----------------------------------------------------------------------------------------------
struct AnaParameters
{
TString outpath;
TString outrootfile;
int events;
int start;
int verbose;
int index{-1};
int pt1{0};
int pt2{0};
int pt3{0};
};
struct tm decodeHadesTimeAndDate(UInt_t date, UInt_t time, bool verbose = false)
{
time_t rawtime;
std::time(&rawtime);
struct tm t = *localtime(&rawtime);
t.tm_mday = (date & 0xFF);
t.tm_mon = 1 + ((date >> 8) & 0xFF) - 1;
t.tm_year = ((date >> 16) & 0xFF);
t.tm_hour = ((time >> 16) & 0xFF) + 1;
t.tm_min = (time >> 8) & 0xFF;
t.tm_sec = time & 0xFF;
if (verbose)
{
printf("Date: %#x Time: %#x\n", date, time);
printf("Decoded %4d-%02d-%02d %02d:%02d:%02d\n", 1900 + t.tm_year, t.tm_mon + 1, t.tm_mday,
t.tm_hour, t.tm_min, t.tm_sec);
}
return t;
}
//##############################################################################################
// sts_drift_time_trend function
//##############################################################################################
Int_t sts_drift_time_trend(HLoop* loop, const AnaParameters& anapars)
{
// Check if loop was properly initialized
// if (!loop->setInput(""))
if (!loop->setInput("-*,+HStsRaw,+HStsCal")) //,+HStart2Hit,+HRpcCal,+HTofRaw"))
//// ,+HFRpcHit,+HFRpcClus
{ // reading file structure
std::cerr << "READBACK: ERROR : cannot read input !\n";
std::exit(EXIT_FAILURE);
}
//----------------------------------------------------------------------------------------------
// Input parameters file
//----------------------------------------------------------------------------------------------
cout << "Input parameters file ...\n";
Int_t mdcMods[6][4] = {{1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1},
{1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}};
HDst::setupSpectrometer("feb22", mdcMods, "sts,frpc");
// pFRpcCalPar->printParams();
//----------------------------------------------------------------------------------------------
// Output parameters file
//----------------------------------------------------------------------------------------------
// Forward Detector Raw data
HCategory* fStsRaw = HCategoryManager::getCategory(catStsRaw, kTRUE, "catStsRaw");
if (!fStsRaw)
{
cout << "No fStsRaw!" << endl;
abort();
}
//----------------------------------------------------------------------------------------------
// Setting parameters for loop over events
//----------------------------------------------------------------------------------------------
Int_t entries = loop->getEntries(); // Number of entries in loop
if (!entries) return -1;
int limit_sta = anapars.start; // Limit START - Where to start the loop
int limit_sto = 0; // Limit STOP - Where to stop the loop
if (anapars.events >= 0)
limit_sto = limit_sta + anapars.events;
else
limit_sto = entries;
if (limit_sto > entries) limit_sto = entries;
//------------------------------------------------------------------------------------------------------
// Important variables definition
//----------------------------------------------------------------------------------------------
char buff1[1000];
char buff2[1000];
Int_t particle_cand_cnt = 0;
Int_t forward_cand_cnt = 0;
Int_t t_PT1, t_PT2, t_PT3;
const auto n_mod = 2;
const auto n_lay = 4;
const auto n_str = 512;
int active[n_mod][n_lay][n_str] = {0};
constexpr int range[n_mod] = {160, 224};
constexpr int urange[n_mod] = {72, 96};
constexpr int ulen[n_mod] = {16, 32};
TFile* outfile = nullptr;
if (anapars.outrootfile.Length()) { outfile = TFile::Open(anapars.outrootfile, "RECREATE"); }
TH1I* h_straws[n_mod][n_lay][n_str];
for (auto i = 0; i < n_mod; ++i)
{
for (auto j = 0; j < n_lay; ++j)
{
for (auto k = 0; k < n_str; ++k)
{
sprintf(buff1, "h_straw_%d_%d_%d", i, j, k);
h_straws[i][j][k] = new TH1I(buff1, buff1, 40, 0, 1000);
}
for (auto k = 0; k < range[i]; ++k)
{
active[i][j][HStsCal::calcStrawIndex(k, 0x03)] = 1;
}
for (auto k = urange[i]; k < urange[i] + ulen[i]; ++k)
{
active[i][j][HStsCal::calcStrawIndex(k, 0x02)] = 1;
}
}
}
//----------------------------------------------------------------------------------------------
// Loop over events in input file
//----------------------------------------------------------------------------------------------
std::string timestamp;
regex str_expr(".*/?[a-z]{2}([0-9]{11})[0-9]{2}.+");
printf("Total events = %d Start = %d To analyze = %d\n", entries, anapars.start,
anapars.events);
//----------------------------------------------------------------------------------------------
// Event loop
for (Int_t i = limit_sta; i < limit_sto; i++)
{
if (i % 10000 == 0)
{
// printf("Event nr.: %d, progress: %.2f%%, time: %f s\n", i, (double)(i - limit_sta) /
// (limit_sto - limit_sta) * 100., timer.RealTime());
printf("Event nr.: %d, Progress: %.2f%%\n", i,
(double)(i - limit_sta) / (limit_sto - limit_sta) * 100.);
}
//----------------------------------------------------------------------------------------------
// Get next event. Categories will be cleared before
loop->nextEvent(i);
HEventHeader* event_header = NULL;
if (!(event_header = gHades->getCurrentEvent()->getHeader())) continue;
// Timestamp pierwszego pliku:
TString new_file;
auto r = loop->isNewFile(new_file);
if (r)
{
auto last = new_file.Last('/');
if (last != kNPOS) {}
std::smatch capture;
string nf = new_file.Data();
if (regex_match(nf, capture, str_expr))
{
char buf[255];
struct tm tm;
memset(&tm, 0, sizeof(tm));
// printf("Captured date: %s\n", capture[1].str().c_str());
strptime(capture[1].str().c_str(), "%y%j%H%M%S", &tm);
strftime(buf, sizeof(buf), "%s", &tm);
timestamp = buf;
}
}
Int_t TBit = (Int_t)event_header->getTBit();
//----------------------------------------------------------------------------------------------
// Checking trigger
bool is_pt1 = false;
bool is_pt2 = false;
bool is_pt3 = false;
// Binary flags for each trigger type
t_PT1 = 0;
t_PT2 = 0;
t_PT3 = 0;
if ((TBit & 2048) == 2048)
{
// bit=1;//PT1
is_pt1 = true;
t_PT1 = 1;
}
if ((TBit & 4096) == 4096)
{
// bit=2;//PT2
is_pt2 = true;
t_PT2 = 1;
}
if ((TBit & 8192) == 8192)
{
// bit=3;//PT3
is_pt3 = true;
t_PT3 = 1;
}
if (fStsRaw)
{
auto raw_cnt = fStsRaw->getEntries();
for (int j = 0; j < raw_cnt; ++j)
{
char m, l, ud;
int s;
HStsRaw* stsraw = HCategoryManager::getObject(stsraw, fStsRaw, j);
stsraw->getAddress(m, l, s, ud);
auto time = stsraw->getTime();
Int_t straw_index = HStsCal::calcStrawIndex(s, ud);
h_straws[m][l][straw_index]->Fill(time);
}
}
} // end eventloop
for (auto i = 0; i < n_mod; ++i)
for (auto j = 0; j < n_lay; ++j)
{
sprintf(buff1, "mod_%d_lay_%d_ts_index%03d.dat", i, j, anapars.index);
std::ofstream ofs(buff1);
ofs << timestamp;
for (auto k = 0; k < n_str; ++k)
{
if (!active[i][j][k]) continue;
ofs << " " << h_straws[i][j][k]->GetBinCenter(h_straws[i][j][k]->GetMaximumBin());
}
ofs << '\n';
}
//------------------------------------------------------------------------------------------------------
// Writing calibration parameters to file and statistics to screen
if (outfile)
{
outfile->Write();
outfile->Close();
}
return 0;
}
int main(int argc, char** argv)
{
TROOT Analysis("Analysis", "compiled analysis macro");
gStyle->SetOptStat(0);
int c;
AnaParameters anapars;
anapars.start = 0;
anapars.events = -1;
anapars.outpath = "";
while (1)
{
static struct option long_options[] = {
/* These options set a flag. */
{"verbose", no_argument, &anapars.verbose, 1},
{"brief", no_argument, &anapars.verbose, 0},
{"pt1", no_argument, &anapars.pt1, 1},
{"pt2", no_argument, &anapars.pt2, 1},
{"pt3", no_argument, &anapars.pt3, 1},
/* These options don’t set a flag.
* We distinguish them by their indices. */
{"dir", required_argument, 0, 'd'},
{"events", required_argument, 0, 'e'},
{"index", required_argument, 0, 'i'},
{"output", required_argument, 0, 'o'},
{"start", required_argument, 0, 's'},
{0, 0, 0, 0}};
/* getopt_long stores the option index here. */
int option_index = 0;
c = getopt_long(argc, argv, "d:e:i:o:s:", long_options, &option_index);
/* Detect the end of the options. */
if (c == -1) break;
switch (c)
{
case 0:
/* If this option set a flag, do nothing else now. */
if (long_options[option_index].flag != 0) break;
printf("option %s", long_options[option_index].name);
if (optarg) printf(" with arg %s", optarg);
printf("\n");
break;
case 'd':
anapars.outpath = optarg;
break;
case 'e':
anapars.events = atol(optarg);
break;
case 'i':
anapars.index = atoi(optarg);
break;
case 'o':
anapars.outrootfile = optarg;
break;
case 's':
anapars.start = atol(optarg);
break;
case '?':
/* getopt_long already printed an error message. */
break;
default:
abort();
}
}
/* Instead of reporting ‘--verbose’
* and ‘--brief’ as they are encountered,
* we report the final status resulting from them. */
if (anapars.verbose) puts("verbose flag is set");
if (anapars.index == -1)
{
std::cerr << "index must be set, use option -i [index] with index >= 0\n";
abort();
}
HLoop* loop = new HLoop(kTRUE);
/* Print any remaining command line arguments (not options). */
if (optind < argc)
{
Bool_t ret;
// printf ("non-option ARGV-elements: ");
while (optind < argc)
{
TString infile = argv[optind++];
if (infile.Contains(","))
ret = loop->addMultFiles(infile);
else if (infile.Contains(".root"))
ret = loop->addFiles(infile);
else
ret = loop->addFilesList(infile);
if (!ret)
{
std::cerr << "READBACK: ERROR : cannot find inputfiles : " << infile.Data() << endl;
std::exit(EXIT_FAILURE);
}
}
}
if (!anapars.pt1 && !anapars.pt2 && !anapars.pt3)
{
anapars.pt1 = 1;
anapars.pt2 = 1;
anapars.pt3 = 1;
}
printf("PT config %d %d %d\n", anapars.pt1, anapars.pt2, anapars.pt3);
if (!anapars.outpath.IsNull()) anapars.outpath += "/";
sts_drift_time_trend(loop, anapars);
exit(0);
}
add_executable(
vertex_monitor
vertex_monitor.cpp
)
target_link_libraries(vertex_monitor FD::Forward ${HYDRA2_LIBRARIES} ${ROOT_LIBRARIES})
# Install the export set for use with the install-tree
install(TARGETS vertex_monitor
RUNTIME
DESTINATION ${CMAKE_INSTALL_BINDIR}
)
install(
FILES
draw_event.C
gui_draw_event.C
DESTINATION
${CMAKE_INSTALL_DATADIR}/forward_tools/macros
)
#include "fwdetdef.h"
#include "hstscal.h"
#include "hcategorymanager.h"
#include "hloop.h"
#include "hruntimedb.h"
#include "hparasciifileio.h"
#include "hstsgeompar.h"
#include "hfrpccalpar.h"
#include "hfrpccluster.h"
#include "hforwardcand.h"
#include <TCanvas.h>
#include <TEllipse.h>
#include <TGraph.h>
#include <TH2.h>
#include <TStyle.h>
#include <TLine.h>
#include <TText.h>
#include <TROOT.h>
void draw_rect(TVirtualPad * pad, float x1, float y1, float x2, float y2) {
pad->cd();
TLine l;
l.SetNDC(0);
l.SetLineStyle(9);
l.DrawLine(x1, y1, x2, y1);
l.DrawLine(x2, y1, x2, y2);
l.DrawLine(x2, y2, x1, y2);
l.DrawLine(x1, y2, x1, y1);
}
void draw_line(TVirtualPad * pad, float x1, float y1, float x2, float y2) {
pad->cd();
TLine l;
l.SetNDC(0);
l.SetLineStyle(9);
l.DrawLine(x1, y1, x2, y2);
}
void draw_straw_circle(TVirtualPad * pad, Float_t x, Float_t z, Int_t color) {
pad->cd();
TEllipse c;
c.SetLineColor(color);
c.DrawEllipse(x, z, 5., 5., 0, 360, 0);
}
void draw_straw_circles(TVirtualPad * pad, TGraph * gr) {
int n = gr->GetN();
double * x = gr->GetX();
double * y = gr->GetY();
for (int i = 0; i < n; ++i)
draw_straw_circle(pad, x[i], y[i], gr->GetMarkerColor());
}
int draw_event(int e, int show_tracks = 1, int select_min_mult = 1) {
gStyle->SetOptStat(0);