Skip to content

Commit

Permalink
remove analysis parameters from filename, and moved them into the roo…
Browse files Browse the repository at this point in the history
…t tree
  • Loading branch information
dorotheapfeiffer committed May 4, 2023
1 parent 5bcdbd0 commit 5f3bbdd
Show file tree
Hide file tree
Showing 6 changed files with 367 additions and 33 deletions.
2 changes: 2 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -88,10 +88,12 @@ target_link_libraries(vmm-sdat ${ROOT_LIBRARIES} fmt pcap)
#---Create a main program using the library
add_executable(convertFile ${CMAKE_CURRENT_SOURCE_DIR}/src/convertFile.cpp)
add_executable(accessTree ${CMAKE_CURRENT_SOURCE_DIR}/analysis/accessTree.cpp)
add_executable(normalizeHisto ${CMAKE_CURRENT_SOURCE_DIR}/analysis/normalizeHisto.cpp)


target_link_libraries(convertFile vmm-sdat)
target_link_libraries(accessTree vmm-sdat)
target_link_libraries(normalizeHisto vmm-sdat)

target_compile_definitions(vmm-sdat PUBLIC __FAVOR_BSD)

Expand Down
6 changes: 6 additions & 0 deletions analysis/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,10 @@ This example can be edited to analyze and plot different data from the tree. In
### Analysis with CINT root script
This example creates a multitude of plots for the LET Multigrid detector that consists of grids and wires. It can be easily adapted to detector with an x/y readout by replacing wires with x and strips with y. To execute the script, the following strips are necessary:
1. start root by typing in a terminal: root
2. in root, type ".x plots_LET.C++("example_LET",".",0)"

### Normalization of histograms
The script normalizeHisto can be used to normalize histograms, or carry out a
This example creates a multitude of plots for the LET Multigrid detector that consists of grids and wires. It can be easily adapted to detector with an x/y readout by replacing wires with x and strips with y. To execute the script, the following strips are necessary:
1. start root by typing in a terminal: root
2. in root, type ".x plots_LET.C++("example_LET",".",0)"
271 changes: 271 additions & 0 deletions analysis/normalizeHisto.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,271 @@
#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include "TFile.h"
#include <vector>
#include <map>
#include "TStyle.h"
#include "TCanvas.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include <TObject.h>
#include <TString.h>
#include "TBufferJSON.h"

#include <chrono>

TString fileInName = "";
TString fileOutName = "";
TString fileNormName = "out";
std::string pAction = "";
std::string pPlot = "";

bool fileInFound = false;
bool fileOutFound = false;
bool fileNormFound = false;

bool plotLog = false;
bool plotPdf = false;
bool plotJSON = false;
int pMap = 56;
int pRebinFactor = 1;
double pZmin = 0;
double pZmax = 0;
double pTimeSignal = 0;
double pTimeBackground = 0;
double pXbins = 640;
double pXmin = 0;
double pXmax=640;
double pYbins = 640;
double pYmin = 0;
double pYmax=640;
std::string pXvar = "pos0";
std::string pYvar = "pos1";
TString cutIn = "";
TString cutNorm = "";

int printUsage(std::string errorMessage, char* argv);

int main(int argc, char**argv) {
std::chrono::time_point<std::chrono::system_clock> timeEnd, timeStart;

if (argc == 1 || argc % 2 == 0) {
return printUsage("Wrong number of arguments!", argv[argc - 1]);
}
for (int i = 1; i < argc; i += 2) {
if (strncmp(argv[i], "-fin", 3) == 0) {
fileInFound = true;
fileInName = argv[i + 1];
} else if (strncmp(argv[i], "-fnorm", 6) == 0) {
fileNormFound = true;
fileNormName = argv[i + 1];
} else if (strncmp(argv[i], "-fout", 5) == 0) {
fileOutFound = true;
fileOutName = argv[i + 1];
} else if (strncmp(argv[i], "-ns", 3) == 0) {
pTimeSignal = atof(argv[i + 1]);
} else if (strncmp(argv[i], "-nb", 3) == 0) {
pTimeBackground = atof(argv[i + 1]);
} else if (strncmp(argv[i], "-act", 4) == 0) {
pAction = argv[i + 1];
} else if (strncmp(argv[i], "-plot", 5) == 0) {
pPlot = argv[i + 1];
if (pPlot.find("j") != std::string::npos) {
plotJSON = true;
}
if (pPlot.find("l") != std::string::npos) {
plotLog = true;
}
if (pPlot.find("p") != std::string::npos) {
plotPdf = true;
}

} else if (strncmp(argv[i], "-map", 4) == 0) {
pMap = atoi(argv[i + 1]);

}
else if (strncmp(argv[i], "-xbin", 5) == 0) {
pXbins = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-xmin", 5) == 0) {
pXmin = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-xmax", 5) == 0) {
pXmax = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-ybin", 5) == 0) {
pYbins = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-ymin", 5) == 0) {
pYmin = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-ymax", 5) == 0) {
pYmax = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-zmin", 5) == 0) {
pZmin = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-zmax", 5) == 0) {
pZmax = atof(argv[i + 1]);
}
else if (strncmp(argv[i], "-xvar", 5) == 0) {
pXvar = argv[i + 1];
}
else if (strncmp(argv[i], "-yvar", 5) == 0) {
pYvar = argv[i + 1];
}
else if (strncmp(argv[i], "-cuti", 5) == 0) {
cutIn = argv[i + 1];
}
else if (strncmp(argv[i], "-cutn", 5) == 0) {
cutNorm = argv[i + 1];
}
else {
return printUsage("Wrong type of argument!", argv[i]);
}
}

if (!fileInFound || !fileInName.EndsWith(".root")) {
return printUsage("Data file has to be loaded with -fin in.root!", nullptr);

}
if (!fileNormFound || !fileNormName.EndsWith(".root")) {
return printUsage("Data file has to be loaded with -fnorm norm.root!", nullptr);

}
if (fileOutName.EndsWith(".root")) {
fileOutName.ReplaceAll(".root", "");
}


timeStart = std::chrono::system_clock::now();
TFile *fIn = new TFile(fileInName, "read");
TH2D *hIn = new TH2D("hIn", "hIn",pXbins,pXmin,pXmax,pYbins, pYmin,pYmax);
TTree *tIn = (TTree*)fIn->Get("clusters_detector");
TString drawChoice = pYvar + ":" + pXvar + ">>hIn";

tIn->Draw(drawChoice, cutIn);


TFile *fNorm = new TFile(fileNormName, "read");
TH2D *hNorm = new TH2D("hNorm", "hNorm",pXbins,pXmin,pXmax,pYbins, pYmin,pYmax);
TTree *tNorm = (TTree*)fNorm->Get("clusters_detector");
drawChoice = pYvar + ":" + pXvar + ">>hNorm";
tNorm->Draw(drawChoice, cutNorm);


TFile *fOut = new TFile(fileOutName+".root", "recreate");
std::cout << fileOutName+".root" << " opened .." << std::endl;

TH2D *hOut = new TH2D("hOut", "hOut",pXbins,pXmin,pXmax,pYbins, pYmin,pYmax);
for (Int_t x = 0; x < pXbins; x++) {
for (Int_t y = 0; y < pYbins; y++) {
if (pTimeSignal * pTimeBackground > 0) {
double valSignal = pTimeBackground * hIn->GetBinContent(x + 1, y + 1);
double valBackground = pTimeSignal * hNorm->GetBinContent(x + 1, y + 1);
if (pAction == "n") {
double val = 0.0;
if (valBackground > 0) {
val = valSignal / valBackground;
}
hOut->SetBinContent(x + 1, y + 1, val);

} else {
double val = 0.0;
if (valSignal - valBackground > 0) {
val = valSignal - valBackground;
}
hOut->SetBinContent(x + 1, y + 1, val);
}
}
else {
double val = hIn->GetBinContent(x + 1, y + 1);
hOut->SetBinContent(x + 1, y + 1, val);
}
}
}
gStyle->SetPalette(pMap);
gStyle->SetPadRightMargin(0.13);

TCanvas *c1 = new TCanvas("c1", "c1", 630, 600);

c1->SetFillColor(0);
c1->SetBorderMode(0);
c1->SetBorderSize(2);
c1->SetFrameBorderMode(0);
hOut->GetXaxis()->SetTitle("x coordinate [mm]");
hOut->GetYaxis()->SetTitle("y coordinate [mm]");
gPad->SetRightMargin(0.18);
hOut->GetXaxis()->SetTitleOffset(1.4);
hOut->GetYaxis()->SetTitleOffset(1.4);
hOut->GetZaxis()->SetTitleOffset(1.5);
hOut->GetZaxis()->SetTitle("counts");
hOut->GetXaxis()->CenterTitle();
hOut->GetYaxis()->CenterTitle();
hOut->GetZaxis()->CenterTitle();
hOut->SetTitle("");
hOut->SetStats(0);
if(pZmax > 0) {
hOut->GetZaxis()->SetRangeUser(pZmin, pZmax);
}
hOut->GetXaxis()->SetRangeUser(pXmin,pXmax);
hOut->GetYaxis()->SetRangeUser(pYmin,pYmax);
hOut->Draw("COLZ");


if (plotJSON) {
TString jsonOut = TBufferJSON::ToJSON(hOut, 3);
std::ofstream fJSON;
fJSON.open(fileOutName+".json", std::ios::out);
fJSON << jsonOut;
fJSON.close();
std::cout << "JSON" << std::endl;
}

if (plotPdf) {
TString pdfName = "";
if (plotLog) {
c1->SetLogz(true);
pdfName = fileOutName+"_log.pdf";
} else {
pdfName = fileOutName+".pdf";
}
c1->SaveAs(pdfName);
}
fOut->cd();
c1->Write();
hOut->Write();
fOut->Close();
fIn->Close();
fNorm->Close();

timeEnd = std::chrono::system_clock::now();

int elapsed_seconds = std::chrono::duration_cast < std::chrono::milliseconds > (timeEnd - timeStart).count();

std::cout << "finished computation in " << elapsed_seconds << " ms\n";

return 0;

}

int printUsage(std::string errorMessage, char* argv) {
if(argv != nullptr) {
std::cout << "\nERROR: " << errorMessage << ": " << argv << std::endl;
}
else {
std::cout << "\nERROR: " << errorMessage << std::endl;
}

printf("\nUsages:\n");
printf("normalize or substract background:\n\t./normalize -fin Run1.root -fnorm Run2.root -fout Out"
" -ns 0 -nb 0 -act normalize\n");

return -1;
}

16 changes: 9 additions & 7 deletions monitoring/monitorData.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
import time as t
import numpy as np
from matplotlib.colors import LogNorm

import gc

#########################
#Edit the paramters below
Expand All @@ -27,7 +27,7 @@
h4_total = np.zeros((channels_y, 128))
h5_total = [0] * channels_x
h6_total = [0] * channels_y
h7_total = [0] * 1000
h7_total = [0] * int(max_charge*0.1)
h8_total = np.zeros((channels_x, channels_y))

h_cluster_rate = [0] * 200
Expand Down Expand Up @@ -85,12 +85,11 @@
print("File " + str(lastFileId) + " already analysed! Waiting for new data..")
t.sleep(2)
continue
args_vmmsdat = ['../build/convertFile', '-f', "./" + name, '-geo', 'example_monitoring.json', '-bc', '40', '-tac', '60', '-th','0', '-cs','1', '-ccs', '3', '-dt', '200', '-mst', '1', '-spc', '500', '-dp', '200', '-coin', 'center-of-masss', '-crl', '0.2', '-cru', '10', '-save', '[[0],[],[0]]', '-algo', '4', '-info', '', '-df','SRS']
args_vmmsdat = ['../build/convertFile', '-f', "./" + name, '-geo', 'example_monitoring.json', '-bc', '40', '-tac', '60', '-th','0', '-cs','1', '-ccs', '3', '-dt', '200', '-mst', '1', '-spc', '500', '-dp', '200', '-coin', 'center-of-masss', '-crl', '0.2', '-cru', '10', '-save', '[[0],[0],[0]]', '-algo', '4', '-info', '', '-df','SRS']
#args_vmmsdat = ['../build/convertFile', '-f', "./" + name, '-vmm', '[[[0,0,0,8],[0,0,0,9],[0,0,0,6],[0,0,0,7],[0,0,0,4],[0,0,0,5],[0,0,0,2],[0,0,0,3],[0,0,0,0],[0,0,0,1],[0,1,1,0],[0,1,1,1],[0,1,1,2],[0,1,1,3],[0,1,1,4],[0,1,1,5],[0,1,1,6],[0,1,1,7],[0,1,1,8],[0,1,1,9]]','-axis', '[[0,0],1],[[0,1],0]', '-bc', '44.02', '-tac', '60', '-th','0', '-cs','1', '-ccs', '2', '-dt', '200', '-mst', '1', '-spc', '500', '-dp', '500', '-coin', 'center-of-masss', '-crl', '0.1', '-cru', '10', '-save', '[[0],[0],[0]]', '-json','0', '-algo', '0', '-info', 'monitor', '-df', 'ESS']

subprocess.call(args_vmmsdat)
lastFileId = fileId
cnt = cnt + 1
now_time = timer()

print("Analysis " + str(cnt) + ": " + str(now_time - start_time) + " s")
Expand Down Expand Up @@ -273,7 +272,7 @@
fig.set_size_inches(fig_w, fig_h)
fig.tight_layout()
fig.savefig("det_stats_" + str(detector_id) + ".png", format="png")

fig.clf()

####################################################################
# Latest histograms
Expand All @@ -296,7 +295,7 @@
h5_total = h5_total + h5[0]
h6 = ax[1, 1].hist(cl['pos1'], bins = channels_y, range = [0.0, channels_y], color=color_y)
h6_total = h6_total + h6[0]
h7 = ax[1, 2].hist(cl['adc0']+cl['adc1'], bins = 1000, range = [0.0, max_charge], color=color_xy)
h7 = ax[1, 2].hist(cl['adc0']+cl['adc1'], bins = int(max_charge*0.1), range = [0.0, max_charge], color=color_xy)
h7_total = h7_total + h7[0]
h8 = ax[1, 3].hist2d(cl['pos0'], cl['pos1'],bins =[channels_x, channels_y],cmap=plt.cm.viridis,range=np.array([(0, channels_x), (0, channels_y)]))
plt.colorbar(h8[3], ax=ax[1, 3], orientation='vertical')
Expand Down Expand Up @@ -334,7 +333,7 @@
fig.set_size_inches(fig_w, fig_h)
fig.tight_layout()
fig.savefig("det_" + str(detector_id) + ".png", format="png")

fig.clf()

####################################################################
# Total histograms
Expand Down Expand Up @@ -395,7 +394,10 @@
####################################################################
# clean up
####################################################################
fig.clf()
plt.close("all")
gc.collect()

end_time = timer()
print("Plot " + str(cnt) + ": " + str(end_time-start_time) + " s")
cnt = cnt + 1
Expand Down
Loading

0 comments on commit 5f3bbdd

Please sign in to comment.