modified: Analyzer.C
modified: GainMatchSX3.C modified: GainMatchSX3Front.C new file: HistPlotter.h
This commit is contained in:
parent
61473ca14e
commit
49de3b64a8
34
Analyzer.C
34
Analyzer.C
|
|
@ -53,8 +53,8 @@ const int MAX_DET = 24;
|
|||
const int MAX_UP = 4;
|
||||
const int MAX_DOWN = 4;
|
||||
const int MAX_BK = 4;
|
||||
double backGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool backGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
double backGain[MAX_DET][MAX_BK] = {{0}};
|
||||
bool backGainValid[MAX_DET][MAX_BK] = {{false}};
|
||||
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
|
||||
|
|
@ -157,35 +157,35 @@ void Analyzer::Begin(TTree * /*tree*/)
|
|||
return;
|
||||
}
|
||||
|
||||
int id, bk, u, d;
|
||||
int id, bk;
|
||||
double gain;
|
||||
while (infile >> id >> bk >> u >> d >> gain)
|
||||
while (infile >> id >> bk >> gain)
|
||||
{
|
||||
backGain[id][bk][u][d] = gain;
|
||||
if (backGain[id][bk][u][d] > 0)
|
||||
backGainValid[id][bk][u][d] = true;
|
||||
backGain[id][bk] = gain;
|
||||
if (backGain[id][bk] > 0)
|
||||
backGainValid[id][bk] = true;
|
||||
else
|
||||
backGainValid[id][bk][u][d] = false;
|
||||
backGainValid[id][bk] = false;
|
||||
}
|
||||
|
||||
infile.close();
|
||||
std::cout << "Loaded back gains from " << filename << std::endl;
|
||||
|
||||
std::string filename = "sx3_GainMatchfront.txt";
|
||||
std::string filename1 = "sx3_GainMatchfront.txt";
|
||||
|
||||
std::ifstream infile(filename);
|
||||
if (!infile.is_open())
|
||||
std::ifstream infile1(filename1);
|
||||
if (!infile1.is_open())
|
||||
{
|
||||
std::cerr << "Error opening " << filename << "!" << std::endl;
|
||||
std::cerr << "Error opening " << filename1 << "!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
int id, bk, u, d;
|
||||
double gain;
|
||||
while (infile >> id >> bk >> u >> d >> gain)
|
||||
int idf, bkf, uf, df;
|
||||
double fgain;
|
||||
while (infile1 >> idf >> bkf >> uf >> df >> fgain)
|
||||
{
|
||||
frontGain[id][bk][u][d] = gain;
|
||||
frontGainValid[id][bk][u][d] = true;
|
||||
frontGain[idf][bkf][uf][df] = fgain;
|
||||
frontGainValid[idf][bkf][uf][df] = true;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -12,6 +12,7 @@
|
|||
#include <algorithm>
|
||||
#include <TProfile.h>
|
||||
#include "Armory/ClassSX3.h"
|
||||
#include "Armory/HistPlotter.h"
|
||||
#include <TGraphErrors.h>
|
||||
#include "TVector3.h"
|
||||
|
||||
|
|
@ -35,13 +36,18 @@ const int MAX_UP = 4;
|
|||
const int MAX_DOWN = 4;
|
||||
const int MAX_BK = 4;
|
||||
|
||||
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
double frontGainUp[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
double frontGainDown[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
TCanvas c("canvas","canvas", 800, 600);
|
||||
|
||||
|
||||
// ==== Configuration Flags ====
|
||||
const bool interactiveMode = false; // If true: show canvas + wait for user
|
||||
const bool interactiveMode = true; // If true: show canvas + wait for user
|
||||
const bool verboseFit = true; // If true: print fit summary and chi²
|
||||
const bool drawCanvases = false; // If false: canvases won't be drawn at all
|
||||
const bool drawCanvases = true; // If false: canvases won't be drawn at all
|
||||
|
||||
// HistPlotter plotter("SX3GainMatchBack.root");
|
||||
|
||||
void GainMatchSX3::Begin(TTree * /*tree*/)
|
||||
{
|
||||
|
|
@ -83,22 +89,26 @@ void GainMatchSX3::Begin(TTree * /*tree*/)
|
|||
}
|
||||
cut1->SetName("UvD");
|
||||
|
||||
// std::string filename = "sx3_GainMatchfront.txt";
|
||||
// plotter.ReadCuts("cuts.txt");
|
||||
|
||||
// std::ifstream infile(filename);
|
||||
// if (!infile.is_open())
|
||||
// {
|
||||
// std::cerr << "Error opening " << filename << "!" << std::endl;
|
||||
// return;
|
||||
// }
|
||||
std::string filename = "sx3_GainMatchfront0.txt";
|
||||
//std::string filename = "sx3_GainMatchfront.txt";
|
||||
|
||||
// int id, bk, u, d;
|
||||
// double gain;
|
||||
// while (infile >> id >> bk >> u >> d >> gain)
|
||||
// {
|
||||
// frontGain[id][bk][u][d] = gain;
|
||||
// frontGainValid[id][bk][u][d] = true;
|
||||
// }
|
||||
std::ifstream infile(filename);
|
||||
if (!infile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening " << filename << "!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
int id, bk, u, d;
|
||||
double gainup,gaindown;
|
||||
while (infile >> id >> bk >> u >> d >> gainup >> gaindown)
|
||||
{
|
||||
frontGainUp[id][bk][u][d] = gainup;
|
||||
frontGainDown[id][bk][u][d] = gaindown;
|
||||
frontGainValid[id][bk][u][d] = true;
|
||||
}
|
||||
}
|
||||
|
||||
Bool_t GainMatchSX3::Process(Long64_t entry)
|
||||
|
|
@ -188,6 +198,7 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
{
|
||||
sx3ChDn = sx3.ch[index];
|
||||
sx3EDn = sx3.e[index];
|
||||
//
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -201,8 +212,9 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
sx3EBk = sx3.e[index];
|
||||
}
|
||||
}
|
||||
sx3EUp*=frontGainUp[sx3ID[i].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2];
|
||||
sx3EDn*=frontGainDown[sx3ID[i].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2];
|
||||
}
|
||||
|
||||
// Only if we found all three channels do we proceed
|
||||
if (sx3ChUp >= 0 && sx3ChDn >= 0 && sx3ChBk >= 0)
|
||||
{
|
||||
|
|
@ -210,18 +222,26 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
hSX3->Fill(sx3ChDn + 4, sx3ChBk);
|
||||
hSX3->Fill(sx3ChUp, sx3ChBk);
|
||||
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
// plotter.Fill2D("hSX3F", 400, 0, 16000, 400, 0, 16000, sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
// Pick detector ID from one of the correlated hits (all same detector)
|
||||
int detID = sx3ID[0].first;
|
||||
|
||||
TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d",
|
||||
detID, sx3ChUp, sx3ChDn, sx3ChBk);
|
||||
TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", detID, sx3ChUp, sx3ChDn, sx3ChBk);
|
||||
TString histName1 = Form("UnCorr_id%d_U%d-D%dvsB%d", detID, sx3ChUp, sx3ChDn, sx3ChBk);
|
||||
|
||||
TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
|
||||
TH2F *hist2d1 = (TH2F *)gDirectory->Get(histName1);
|
||||
if (!hist2d)
|
||||
{
|
||||
hist2d = new TH2F(histName, histName,
|
||||
400, 0, 16000, 400, 0, 16000);
|
||||
}
|
||||
if (!hist2d1)
|
||||
{
|
||||
hist2d1 = new TH2F(histName1, histName1,
|
||||
800, -1, 1, 800, 0, 4000);
|
||||
}
|
||||
|
||||
if (sx3EBk > 100 || sx3EUp > 100 || sx3EDn > 100)
|
||||
{
|
||||
|
|
@ -233,6 +253,8 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
}
|
||||
|
||||
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
hist2d1->Fill((sx3EUp - sx3EDn) / (sx3EUp + sx3EDn), sx3EBk);
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -246,7 +268,7 @@ void GainMatchSX3::Terminate()
|
|||
double backSlope[MAX_DET][MAX_BK] = {{0}};
|
||||
bool backSlopeValid[MAX_DET][MAX_BK] = {{false}};
|
||||
|
||||
std::ofstream outFile("sx3_BackGains.txt");
|
||||
std::ofstream outFile("sx3_BackGains1.txt");
|
||||
if (!outFile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening sx3_BackGains.txt for writing!" << std::endl;
|
||||
|
|
@ -295,7 +317,6 @@ void GainMatchSX3::Terminate()
|
|||
|
||||
if (drawCanvases)
|
||||
{
|
||||
TCanvas *c = new TCanvas(Form("c_%d_%d", id, bk), "Back Fit", 800, 600);
|
||||
g.SetTitle(Form("Detector %d Back %d: (Up+Dn) vs Back", id, bk));
|
||||
g.SetMarkerStyle(20);
|
||||
g.SetMarkerColor(kBlue);
|
||||
|
|
@ -316,12 +337,12 @@ void GainMatchSX3::Terminate()
|
|||
|
||||
if (interactiveMode)
|
||||
{
|
||||
c->Update();
|
||||
c.Update();
|
||||
gPad->WaitPrimitive();
|
||||
}
|
||||
else
|
||||
{
|
||||
c->Close();
|
||||
c.Close();
|
||||
}
|
||||
}
|
||||
else
|
||||
|
|
@ -329,7 +350,7 @@ void GainMatchSX3::Terminate()
|
|||
g.Fit(&f, "QNR");
|
||||
}
|
||||
|
||||
double slope = 1/f.GetParameter(0);
|
||||
double slope = 1 / f.GetParameter(0);
|
||||
if (std::abs(slope - 1.0) < 0.3) // sanity check
|
||||
{
|
||||
backSlope[id][bk] = slope;
|
||||
|
|
@ -353,6 +374,9 @@ void GainMatchSX3::Terminate()
|
|||
600, 0, 16000, 600, 0, 16000);
|
||||
TH2F *hAsym = new TH2F("hAsym", "Up vs Dn divide corrected back;Up/Back E;Dn/Back E",
|
||||
400, 0.0, 1.0, 400, 0.0, 1.0);
|
||||
TH2F *hAsymUnorm = new TH2F("hAsymUnorm", "Up vs Dn;Up E;Dn E",
|
||||
800, 0.0, 4000.0, 800, 0.0, 4000.0);
|
||||
|
||||
|
||||
// Fill histograms using corrected back energies
|
||||
for (const auto &kv : dataPoints)
|
||||
|
|
@ -375,8 +399,20 @@ void GainMatchSX3::Terminate()
|
|||
double correctedBack = eBk * slope;
|
||||
double asym = (eUp - eDn) / updn;
|
||||
|
||||
hFVB->Fill(updn,correctedBack );
|
||||
hFVB->Fill(updn, correctedBack);
|
||||
hAsym->Fill(eUp / correctedBack, eDn / correctedBack);
|
||||
hAsymUnorm->Fill(eUp, eDn);
|
||||
TString histNamex = Form("CorrBack_id%d_U%d-D%dvsB%d", id, u, d, bk);
|
||||
|
||||
TH2F *hist2dx = (TH2F *)gDirectory->Get(histNamex);
|
||||
if (!hist2dx)
|
||||
{
|
||||
hist2dx = new TH2F(histNamex, histNamex,
|
||||
800, -1, 1, 800, 0, 4000);
|
||||
}
|
||||
|
||||
hist2dx->Fill((eUp - eDn) / (eUp + eDn), correctedBack);
|
||||
}
|
||||
}
|
||||
// plotter.FlushToDisk();
|
||||
}
|
||||
|
|
|
|||
|
|
@ -30,6 +30,7 @@ SX3 sx3_contr;
|
|||
TCutG *cut;
|
||||
TCutG *cut1;
|
||||
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
|
||||
TCanvas c(Form("canvas"), "Fit", 800, 600);
|
||||
|
||||
// Gain arrays
|
||||
|
||||
|
|
@ -43,9 +44,9 @@ double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
|||
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
double uvdslope[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
// ==== Configuration Flags ====
|
||||
const bool interactiveMode = false; // If true: show canvas + wait for user
|
||||
const bool interactiveMode = true; // If true: show canvas + wait for user
|
||||
const bool verboseFit = true; // If true: print fit summary and chi²
|
||||
const bool drawCanvases = false; // If false: canvases won't be drawn at all
|
||||
const bool drawCanvases = true; // If false: canvases won't be drawn at all
|
||||
|
||||
void GainMatchSX3Front::Begin(TTree * /*tree*/)
|
||||
{
|
||||
|
|
@ -204,14 +205,14 @@ Bool_t GainMatchSX3Front::Process(Long64_t entry)
|
|||
// Fill the histogram for the front vs back
|
||||
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
if (sx3.e[i] > 100) // && sx3.id[i] == 4)
|
||||
if (sx3.e[i] > 100 && sx3.id[i] == 3)
|
||||
{
|
||||
// back gain correction
|
||||
|
||||
// Fill the histogram for the front vs back with gain correction
|
||||
hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
// Fill the index vs energy histogram
|
||||
hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
|
||||
// hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
// // Fill the index vs energy histogram
|
||||
// hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
|
||||
// }
|
||||
// {
|
||||
TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk);
|
||||
|
|
@ -223,8 +224,8 @@ Bool_t GainMatchSX3Front::Process(Long64_t entry)
|
|||
|
||||
// if (sx3ChBk == 2)
|
||||
// printf("Found back channel Det %d Back %d \n", sx3.id[i], sx3ChBk);
|
||||
// hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
|
||||
// hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
|
||||
hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
|
|
@ -292,7 +293,7 @@ void GainMatchSX3Front::Terminate()
|
|||
// f.SetParameter(0, 1.0); // Initial guess for the gain
|
||||
// g.Fit(&f, "R");
|
||||
|
||||
const double fixedError = 20.0; // in ADC channels
|
||||
const double fixedError = 0.0; // in ADC channels
|
||||
|
||||
std::vector<double> xVals, yVals, exVals, eyVals;
|
||||
|
||||
|
|
@ -305,18 +306,17 @@ void GainMatchSX3Front::Terminate()
|
|||
xVals.push_back(x);
|
||||
yVals.push_back(y);
|
||||
exVals.push_back(fixedError); // error in up energy
|
||||
// eyVals.push_back(fixedError); // error in down energy
|
||||
eyVals.push_back(0.); // error in down energy
|
||||
}
|
||||
|
||||
// Build TGraphErrors with errors
|
||||
TGraphErrors g(xVals.size(), xVals.data(), yVals.data(), exVals.data(), eyVals.data());
|
||||
|
||||
TF1 f("f", "[0]*x", 0, 16000);
|
||||
TF1 f("f", "[0]*x+[1]", 0, 16000);
|
||||
f.SetParameter(0, -1.0); // Initial guess
|
||||
|
||||
if (drawCanvases)
|
||||
{
|
||||
TCanvas *c = new TCanvas(Form("c_%d_%d_%d_%d", id, bk, u, d), "Fit", 800, 600);
|
||||
g.SetTitle(Form("Detector %d: U%d D%d B%d", id, u, d, bk));
|
||||
g.SetMarkerStyle(20);
|
||||
g.SetMarkerColor(kBlue);
|
||||
|
|
@ -337,12 +337,12 @@ void GainMatchSX3Front::Terminate()
|
|||
|
||||
if (interactiveMode)
|
||||
{
|
||||
c->Update();
|
||||
c.Update();
|
||||
gPad->WaitPrimitive();
|
||||
}
|
||||
else
|
||||
{
|
||||
c->Close(); // Optionally avoid clutter in batch
|
||||
c.Close(); // Optionally avoid clutter in batch
|
||||
}
|
||||
}
|
||||
else
|
||||
|
|
@ -350,11 +350,23 @@ void GainMatchSX3Front::Terminate()
|
|||
g.Fit(&f, "QNR");
|
||||
}
|
||||
|
||||
frontGain[id][bk][u][d] = f.GetParameter(0);
|
||||
frontGainValid[id][bk][u][d] = true;
|
||||
double slope = f.GetParameter(0);
|
||||
double intercept = f.GetParameter(1);
|
||||
|
||||
outFile << id << " " << bk << " " << u << " " << d << " " << frontGain[id][bk][u][d] << std::endl;
|
||||
// printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]);
|
||||
if (std::abs(slope + 1.0) < 0.3) // sanity check
|
||||
{
|
||||
frontGain[id][bk][u][d] = slope;
|
||||
|
||||
frontGainValid[id][bk][u][d] = true;
|
||||
outFile << id << " " << bk << " " << u << " " << d << " " << TMath::Abs(slope)/intercept << " " << 1.0/intercept << std::endl;
|
||||
printf("Back slope Det%d Bk%d → %.4f\n", id, bk, slope);
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cerr << "Warning: Bad slope for Det" << id << " Bk" << bk
|
||||
<< " slope=" << f.GetParameter(0) << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
outFile.close();
|
||||
|
|
@ -369,12 +381,10 @@ void GainMatchSX3Front::Terminate()
|
|||
|
||||
auto [id, bk, u, d] = kv.first;
|
||||
double front;
|
||||
if (abs((uvdslope[id][bk][u][d] + 1) < 0.2))
|
||||
{
|
||||
if (frontGainValid[id][bk][u][d])
|
||||
front = frontGain[id][bk][u][d];
|
||||
}
|
||||
else
|
||||
front=0.;
|
||||
continue;
|
||||
for (const auto &pr : kv.second)
|
||||
{
|
||||
double eBk, eUp, eDn;
|
||||
|
|
|
|||
414
HistPlotter.h
Normal file
414
HistPlotter.h
Normal file
|
|
@ -0,0 +1,414 @@
|
|||
#ifndef HISTPLOTTER_H
|
||||
#define HISTPLOTTER_H
|
||||
#include <TCanvas.h>
|
||||
#include <TROOT.h>
|
||||
#include <TSystem.h>
|
||||
#include <TStyle.h>
|
||||
#include <iostream>
|
||||
#include <TFile.h>
|
||||
#include <TMemFile.h>
|
||||
#include <TH1.h>
|
||||
#include <TH2.h>
|
||||
#include <TCutG.h>
|
||||
#include <signal.h>
|
||||
#include <cstdlib>
|
||||
#include <utility>
|
||||
#include <fstream>
|
||||
#include <sstream>
|
||||
#include <unordered_map>
|
||||
#include <set>
|
||||
#include <TGraphErrors.h>
|
||||
|
||||
class HistPlotter {
|
||||
private:
|
||||
long long barrier_count, barrier_limit; //meant to keep track of how often to call FillN() on histograms
|
||||
enum {TFILE, TMEMFILE} filetype;
|
||||
std::unordered_map<std::string,TObject*> oMap; //!< Maps std::string to all TH1, TH2 objects in the class
|
||||
std::unordered_map<std::string,TObject*> cutsMap; //!< Maps std::string to TCutG objects held by the class
|
||||
std::set<std::string> folderList; //!< List of all folder names used to nest objects
|
||||
std::unordered_map<TObject*,std::string> foldersForObjects; //!< Map that returns the folder corresponding to the object whose pointer is specified
|
||||
TFile *ofile=nullptr; //!< TFile pointer for the output file
|
||||
TMemFile *omfile=nullptr; //!< TFile pointer for the output memfile
|
||||
|
||||
//Caches to permit FillN() calls
|
||||
std::unordered_map<std::string, std::vector<double>> onedimcache;
|
||||
std::unordered_map<std::string, std::pair<std::vector<double>, std::vector<double>>> twodimcache;
|
||||
inline void FillN_All_Histograms();
|
||||
public:
|
||||
HistPlotter(std::string outfile, std::string type);
|
||||
inline void FlushToDisk(); //!< Writes all objects to file before closing, nesting objects in folders as is found necessary
|
||||
inline void PrintObjects(); //!< Dump objects to std::cout for inspection
|
||||
inline void ReadCuts(std::string);
|
||||
inline TCutG* FindCut(std::string cut) {
|
||||
return static_cast<TCutG*>(cutsMap.at(cut));
|
||||
}
|
||||
inline void set_barrier_limit(long long limit) { barrier_limit = limit; }
|
||||
inline void barrier_increment() {
|
||||
barrier_count++;
|
||||
if(barrier_count == barrier_limit) {
|
||||
FillN_All_Histograms();
|
||||
barrier_count=0;
|
||||
}
|
||||
}
|
||||
/*! \fn void FindCut()
|
||||
\brief
|
||||
- Searches for a cut by name 'cut' in the internal list of cuts 'cutsMap'. Ugly fails (via unresolved at()) if such a cut isn't found.
|
||||
\param filename - name of the plainxtext file containing the cut file locations and identifiers
|
||||
\return Pointer to the TCutG object that matches the name. Very useful to use this as plotter.FindCut("protonbarrelpid")->IsInside(deltaE, E) for instance.
|
||||
*/
|
||||
|
||||
inline void SetNewTitle(std::string name, std::string title) {
|
||||
auto result = oMap.find(name); //result is an iterator
|
||||
if(result==oMap.end()) return; //no warnings, could be changed in future
|
||||
else
|
||||
static_cast<TNamed*>(oMap.at(name))->SetTitle(title.c_str()); // set new title
|
||||
}
|
||||
|
||||
//Smart functions that create a new histogram if it doesn't exist.
|
||||
inline void FillGraph(const std::string &name, float valuex, float valuey, float errx=0, float erry=0);
|
||||
inline void Fill1D(const std::string& name,int nbinsx, float xlow, float xhigh, float value);
|
||||
inline void Fill2D(const std::string& name,int nbinsx, float xlow, float xhigh
|
||||
,int nbinsy, float ylow, float yhigh, float valuex, float valuey);
|
||||
inline void Fill1D(const std::string& name,int nbinsx, float xlow, float xhigh, float value, const std::string& folder);
|
||||
inline void Fill2D(const std::string& name,int nbinsx, float xlow, float xhigh
|
||||
,int nbinsy, float ylow, float yhigh, float valuex, float valuey, const std::string& folder);
|
||||
//TObject* findObject(std::string key);
|
||||
};
|
||||
|
||||
HistPlotter::HistPlotter(std::string outfile, std::string type="") {
|
||||
/*!
|
||||
\brief Constructor. Opens a TFile instance with the specified filename
|
||||
\param outfile : std::string that holds the desired output ROOT filename
|
||||
\return None
|
||||
*/
|
||||
if(type=="" || type == "TFILE") {
|
||||
ofile = new TFile(outfile.c_str(),"recreate");
|
||||
filetype = TFILE;
|
||||
} else if(type =="TMEMFILE") {
|
||||
omfile = new TMemFile(outfile.c_str(),"recreate");
|
||||
filetype=TMEMFILE;
|
||||
} else {
|
||||
std::cout << "Unknown type "<< type << " specified for HistPlotter (use \"TFILE\" or \"TMEMFILE\"), using default \"TFILE\" " << std::endl;
|
||||
ofile = new TFile(outfile.c_str(),"recreate");
|
||||
filetype = TFILE;
|
||||
}
|
||||
barrier_count=0;
|
||||
barrier_limit=1000;
|
||||
}
|
||||
|
||||
void HistPlotter::FillN_All_Histograms() {
|
||||
for(auto it=oMap.begin(); it!=oMap.end(); it++ ) {
|
||||
//it->first is std::string 'name', it->second is the TObject
|
||||
if(it->second->InheritsFrom("TH1F")) {
|
||||
//FillN(size, array-of-doubles, array-of-weights); //we set array-of-weights to (1,1,1,.. (size)
|
||||
static_cast<TH1F*>(it->second)->FillN(onedimcache[it->first].size(), //size
|
||||
onedimcache[it->first].data(), //array
|
||||
std::vector<double>(onedimcache[it->first].size(),1.0).data()); //weight of ones
|
||||
onedimcache[it->first].clear();
|
||||
} else if(it->second->InheritsFrom("TH2F")) {
|
||||
//FillN(size, array-of-doubles, array-of-weights); //we set array-of-weights to (1,1,1,.. (size))
|
||||
static_cast<TH2F*>(it->second)->FillN(twodimcache[it->first].first.size(), //size
|
||||
twodimcache[it->first].first.data(), //x array
|
||||
twodimcache[it->first].second.data(), //y array
|
||||
std::vector<double>(twodimcache[it->first].first.size(),1.0).data()); //weight of ones
|
||||
twodimcache[it->first].first.clear();
|
||||
twodimcache[it->first].second.clear();
|
||||
}
|
||||
}
|
||||
std::cout << "." << std::endl;
|
||||
}
|
||||
void HistPlotter::FlushToDisk() {
|
||||
/*! \fn void FlushToDisk()
|
||||
\brief Function that can be used at any point to exit smoothly by saving all ROOT objects in memory
|
||||
to the output file before closing it. Obeys the binding of histograms to separate folders, if so specified.
|
||||
\return No return -- void
|
||||
*/
|
||||
if(filetype==TMEMFILE && omfile) {
|
||||
std::cout << "Not flushing a TMemfile .. exiting .." << std::endl;
|
||||
delete omfile;
|
||||
return;
|
||||
}
|
||||
if(ofile->IsZombie() || !ofile) {
|
||||
std::cerr << "Output file is zombie, finishing up without writing to disk!" << std::endl;
|
||||
return;
|
||||
}
|
||||
FillN_All_Histograms();
|
||||
for(auto it=oMap.begin(); it!=oMap.end(); it++ ) {
|
||||
//omap maps: name(first) to object address(second).
|
||||
// foldersForObjects maps: object address(first) to foldername(second)
|
||||
auto result = foldersForObjects.find(it->second); //returns <TObject* histogram,std::string foldername> pair if found
|
||||
if(result!=foldersForObjects.end()) { //we try to create folder if needed and cd to it
|
||||
ofile->mkdir(result->second.c_str(),"",kTRUE); // args: name, title, returnExistingDirectory
|
||||
ofile->cd(result->second.c_str());
|
||||
} else {
|
||||
ofile->cd(); //toplevel for all default histograms. Default setting
|
||||
}
|
||||
it->second->Write();
|
||||
}
|
||||
|
||||
//Create a directory for all cuts, and save all cuts in them
|
||||
ofile->mkdir("gCUTS","",kTRUE);
|
||||
ofile->cd("gCUTS");
|
||||
for(auto it=cutsMap.begin(); it!=cutsMap.end(); it++) {
|
||||
(static_cast<TNamed*>(it->second))->SetName(it->first.c_str());
|
||||
it->second->Write();
|
||||
}
|
||||
ofile->Close();
|
||||
std::cout << "Wrote " << oMap.size() << " histograms to TFile " << std::string(ofile->GetName()) << std::endl;
|
||||
}
|
||||
|
||||
void HistPlotter::FillGraph(const std::string& name, float valuex, float valuey, float errx, float erry) {
|
||||
/*! \fn void FillGraph()
|
||||
\brief
|
||||
- Creates a TGraphError in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
|
||||
- Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TGraph*
|
||||
|
||||
\param name name of the TGraph
|
||||
\param valuex The xvalue
|
||||
\param valuey The yvalue
|
||||
\param errx The x error
|
||||
\param erry The y error
|
||||
\return No return void
|
||||
*/
|
||||
auto result = oMap.find(name);
|
||||
if(result==oMap.end()) {
|
||||
TGraphErrors *tempG = new TGraphErrors();
|
||||
tempG->SetName(name.c_str());
|
||||
oMap.insert(std::make_pair(name,static_cast<TObject*>(tempG)));
|
||||
}
|
||||
if(!oMap.at(name)->InheritsFrom("TGraphErrors")) {
|
||||
std::cerr << "Object " << name << " refers to something other than a TGraph*, not filling it hence!" << std::endl;
|
||||
std::cerr << "Abort.." << std::endl;
|
||||
FlushToDisk();
|
||||
exit(-1);
|
||||
}
|
||||
// static_cast<TGraphErrors*>(oMap.at(name))->AddPointError(valuex,valuey,errx,erry);
|
||||
}
|
||||
|
||||
void HistPlotter::Fill1D(const std::string& name, int nbinsx, float xlow, float xhigh, float value) {
|
||||
/*! \fn void Fill1D()
|
||||
\brief
|
||||
- Creates a TH1F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
|
||||
- Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH1*
|
||||
|
||||
\param name name of the TH1F histogram
|
||||
\param nbinsx Number of bins in the histogram
|
||||
\param xlow Lower limit on x-axis
|
||||
\param xhigh Upper limit on x-axis
|
||||
\param value The bin corresponding to value in (nbinsx, xlow, xhigh) is incremented by 1
|
||||
\return No return void
|
||||
*/
|
||||
auto result = oMap.find(name); //result is an iterator
|
||||
if(result==oMap.end()) {
|
||||
TH1F* temp1D = new TH1F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh);
|
||||
oMap.insert(std::make_pair(name,static_cast<TObject*>(temp1D)));
|
||||
onedimcache.insert(std::make_pair(name, std::vector<double>()));
|
||||
onedimcache[name].reserve(16384);
|
||||
} else if(foldersForObjects.find(oMap.at(name))!=foldersForObjects.end()) { //shouldn't have a folder associated with it
|
||||
std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] << ", choose a different name for the histogram to be stored in toplevel .." << std::endl;
|
||||
}
|
||||
|
||||
//Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
|
||||
if(!oMap.at(name)->InheritsFrom("TH1F")) {
|
||||
std::cerr << "Object " << name << " refers to something other than a TH1*, not filling it hence!" << std::endl;
|
||||
std::cerr << "Abort.." << std::endl;
|
||||
FlushToDisk();
|
||||
exit(-1);
|
||||
}
|
||||
onedimcache[name].emplace_back(value);
|
||||
//static_cast<TH1F*>(oMap.at(name))->Fill(value);
|
||||
}
|
||||
|
||||
void HistPlotter::Fill1D(const std::string& name, int nbinsx, float xlow, float xhigh, float value, const std::string& foldername) {
|
||||
/*! \fn void Fill1D()
|
||||
\brief
|
||||
- Creates a TH1F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
|
||||
- Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH1*
|
||||
- Remembers the foldername this particular histogram maps to, if provided. If not, defaults to toplevel.
|
||||
|
||||
\param name name of the TH1F histogram
|
||||
\param nbinsx Number of bins in the histogram
|
||||
\param xlow Lower limit on x-axis
|
||||
\param xhigh Upper limit on x-axis
|
||||
\param value The bin corresponding to value in (nbinsx, xlow, xhigh) is incremented by 1
|
||||
\param foldername Name of the folder to put this histogram into. Defaults to toplevel if left empty
|
||||
\return No return -- void
|
||||
*/
|
||||
|
||||
auto result = oMap.find(name); //result is an iterator
|
||||
if(result==oMap.end()) {
|
||||
TH1F* temp1D = new TH1F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh);
|
||||
oMap.insert(std::make_pair(name,static_cast<TObject*>(temp1D)));
|
||||
onedimcache.insert(std::make_pair(name, std::vector<double>()));
|
||||
onedimcache[name].reserve(16384);
|
||||
if(foldername!="") {
|
||||
if(folderList.find(foldername)==folderList.end()) {
|
||||
folderList.insert(foldername);
|
||||
}
|
||||
foldersForObjects.insert(std::make_pair(static_cast<TObject*>(temp1D),foldername));
|
||||
}
|
||||
} else {
|
||||
//object is present in map, but we enforce unique names
|
||||
//it must already have a folder attached to it
|
||||
if(foldersForObjects.find(oMap.at(name))==foldersForObjects.end()) {
|
||||
std::cerr << "Object " << name << " already registered at toplevel, choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
|
||||
} else if(foldersForObjects[oMap[name]]!=foldername) {
|
||||
std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] << ", choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
|
||||
}
|
||||
}
|
||||
//Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
|
||||
if(!oMap.at(name)->InheritsFrom("TH1F")) {
|
||||
std::cerr << "Object " << name << " refers to something other than a TH1*, not filling it hence!" << std::endl;
|
||||
std::cerr << "Abort.." << std::endl;
|
||||
FlushToDisk();
|
||||
exit(-1);
|
||||
}
|
||||
onedimcache[name].emplace_back(value);
|
||||
//static_cast<TH1F*>(oMap.at(name))->Fill(value);
|
||||
}
|
||||
|
||||
void HistPlotter::Fill2D(const std::string& name, int nbinsx, float xlow, float xhigh, int nbinsy, float ylow, float yhigh, float valuex, float valuey) {
|
||||
/*! \fn void Fill2D()
|
||||
\brief
|
||||
- Creates a TH2F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
|
||||
- Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH2*
|
||||
\param name name of the TH1F histogram
|
||||
\param nbinsx Number of xbins in the histogram
|
||||
\param xlow Lower limit on x-axis
|
||||
\param xhigh Upper limit on x-axis
|
||||
\param nbinsy Number of ybins in the histogram
|
||||
\param ylow Lower limit on y-axis
|
||||
\param yhigh Upper limit on y-axis
|
||||
\param valuex
|
||||
\param valuey The bin corresponding to (valuex, valuey) in (nbinsx, xlow, xhigh, ybinsx, ylow, yhigh) is incremented by 1
|
||||
\return No return -- void
|
||||
*/
|
||||
|
||||
auto result = oMap.find(name); //result is an iterator
|
||||
if(result==oMap.end()) {
|
||||
TH2F* temp2D = new TH2F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh, nbinsy, ylow, yhigh);
|
||||
oMap.insert(std::make_pair(name,static_cast<TObject*>(temp2D)));
|
||||
twodimcache.insert(std::make_pair(name, std::make_pair(std::vector<double>(),std::vector<double>())));
|
||||
twodimcache[name].first.reserve(16384);
|
||||
twodimcache[name].second.reserve(16384);
|
||||
} else if(foldersForObjects.find(oMap.at(name))!=foldersForObjects.end()) { //shouldn't have a folder associated with it
|
||||
std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] << ", choose a different name for the histogram to be stored in toplevel .." << std::endl;
|
||||
}
|
||||
|
||||
//Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
|
||||
if(!oMap.at(name)->InheritsFrom("TH2F")) {
|
||||
std::cerr << "Object " << name << " refers to something other than a TH2*, not filling it hence!" << std::endl;
|
||||
std::cerr << "Abort.." << std::endl;
|
||||
FlushToDisk();
|
||||
exit(-1);
|
||||
}
|
||||
twodimcache[name].first.emplace_back(valuex);
|
||||
twodimcache[name].second.emplace_back(valuey);
|
||||
//static_cast<TH2F*>(oMap.at(name))->Fill(valuex,valuey);
|
||||
}
|
||||
|
||||
void HistPlotter::Fill2D(const std::string& name, int nbinsx, float xlow, float xhigh, int nbinsy, float ylow, float yhigh, float valuex, float valuey, const std::string& foldername) {
|
||||
/*! \fn void Fill2D()
|
||||
\brief
|
||||
- Creates a TH2F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
|
||||
- Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH2*
|
||||
- Remembers the foldername this particular histogram maps to, if provided. If not defaults to toplevel
|
||||
|
||||
\param name name of the TH1F histogram
|
||||
\param nbinsx Number of xbins in the histogram
|
||||
\param xlow Lower limit on x-axis
|
||||
\param xhigh Upper limit on x-axis
|
||||
\param nbinsy Number of ybins in the histogram
|
||||
\param ylow Lower limit on y-axis
|
||||
\param yhigh Upper limit on y-axis
|
||||
\param valuex
|
||||
\param valuey The bin corresponding to (valuex, valuey) in (nbinsx, xlow, xhigh, ybinsx, ylow, yhigh) is incremented by 1
|
||||
\param foldername Name of the folder to put this histogram into. Defaults to toplevel if left empty
|
||||
\return No return -- void
|
||||
*/
|
||||
|
||||
auto result = oMap.find(name); //result is an iterator
|
||||
if(result==oMap.end()) {
|
||||
TH2F* temp2D = new TH2F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh, nbinsy, ylow, yhigh);
|
||||
oMap.insert(std::make_pair(name,static_cast<TObject*>(temp2D)));
|
||||
twodimcache.insert(std::make_pair(name, std::make_pair(std::vector<double>(),std::vector<double>())));
|
||||
twodimcache[name].first.reserve(16384);
|
||||
twodimcache[name].second.reserve(16384);
|
||||
if(foldername!="") {
|
||||
if(folderList.find(foldername)==folderList.end()) {
|
||||
folderList.insert(foldername);
|
||||
}
|
||||
foldersForObjects.insert(std::make_pair(static_cast<TObject*>(temp2D),foldername));
|
||||
}
|
||||
} else {
|
||||
//object is present in map, but we enforce unique names
|
||||
//it must already have a folder attached to it
|
||||
if(foldersForObjects.find(oMap.at(name))==foldersForObjects.end()) {
|
||||
std::cerr << "Object " << name << " already registered at toplevel, choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
|
||||
} else if(foldersForObjects[oMap.at(name)]!=foldername) {
|
||||
std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] << ", choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
//Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
|
||||
if(!oMap.at(name)->InheritsFrom("TH2F")) {
|
||||
std::cerr << "Object " << name << " refers to something other than a TH2*, not filling it hence!" << std::endl;
|
||||
std::cerr << "Abort.." << std::endl;
|
||||
FlushToDisk();
|
||||
exit(-1);
|
||||
}
|
||||
twodimcache[name].first.emplace_back(valuex);
|
||||
twodimcache[name].second.emplace_back(valuey);
|
||||
//static_cast<TH2F*>(oMap.at(name))->Fill(valuex,valuey);
|
||||
}
|
||||
|
||||
void HistPlotter::ReadCuts(std::string filename) {
|
||||
/*! \fn void ReadCuts()
|
||||
\brief Reads a list of cuts from a file. The file must have the format below, two columns
|
||||
- Column#1 - path to a file that contains a single TCutG object named "CUTG", the default name in ROOT.
|
||||
- Column#2 - The identifier name you plan to use in the code, like 'protonbarrelpid' or something, that will be searched by FindCut()
|
||||
\param filename name of the plainxtext file containing the cut file locations and identifiers
|
||||
\return No return -- void
|
||||
*/
|
||||
|
||||
std::ifstream infile;
|
||||
infile.open(filename);
|
||||
std::string cutfilename, cutname;
|
||||
for(std::string line; std::getline(infile, line); ) {
|
||||
if(line.size()!=0 && line[0]=='#')
|
||||
; //don't do anything with '#' lines
|
||||
else {
|
||||
std::stringstream ss(line);
|
||||
ss>>cutfilename>>cutname;
|
||||
|
||||
TFile f(cutfilename.c_str());
|
||||
if(f.IsZombie()) {
|
||||
std::cerr << "Cannot open cutfile " << cutfilename << " .. skipping.." << std::endl;
|
||||
continue;
|
||||
}
|
||||
TCutG *cut = (TCutG*)(f.Get("CUTG"));
|
||||
cutsMap.insert(std::make_pair(cutname,static_cast<TObject*>(cut)));
|
||||
f.Close();
|
||||
} //else
|
||||
}//for loop
|
||||
infile.close();
|
||||
}
|
||||
|
||||
void HistPlotter::PrintObjects() {
|
||||
/*
|
||||
void PrintObjects()
|
||||
Prints the contents of the unordered_maps oMap and cutsMap to facilitate debugging
|
||||
|
||||
*/
|
||||
std::cout << "Type | Name " << std::endl;
|
||||
std::cout << "---- | --------------------- " << std::endl;
|
||||
for(auto it=oMap.begin(); it!=oMap.end(); it++ ) {
|
||||
std::cout << it->second->ClassName() << " | "<< it->first << std::endl;
|
||||
}
|
||||
for(auto it=cutsMap.begin(); it!=cutsMap.end(); it++ ) {
|
||||
std::cout << it->second->ClassName() << " | "<< it->first << std::endl;
|
||||
}
|
||||
std::cout << "---- | --------------------- " << std::endl;
|
||||
}
|
||||
|
||||
#endif
|
||||
Loading…
Reference in New Issue
Block a user