diff --git a/Analyzer.C b/Analyzer.C index 2c1ea08..a2eb4d3 100644 --- a/Analyzer.C +++ b/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; } } diff --git a/GainMatchSX3.C b/GainMatchSX3.C index 3f58e2c..3e40eec 100644 --- a/GainMatchSX3.C +++ b/GainMatchSX3.C @@ -12,6 +12,7 @@ #include #include #include "Armory/ClassSX3.h" +#include "Armory/HistPlotter.h" #include #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 verboseFit = true; // If true: print fit summary and chi² -const bool drawCanvases = false; // If false: canvases won't be drawn at all +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 = 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; @@ -283,7 +305,7 @@ void GainMatchSX3::Terminate() continue; // not enough statistics // Build graph with errors - const double fixedError = 0.0; // ADC channels + const double fixedError = 0.0; // ADC channels std::vector exVals(udE.size(), 0.0); // no x error std::vector eyVals(udE.size(), fixedError); // constant y error @@ -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(); } diff --git a/GainMatchSX3Front.C b/GainMatchSX3Front.C index 61de72e..9166c2c 100644 --- a/GainMatchSX3Front.C +++ b/GainMatchSX3Front.C @@ -30,6 +30,7 @@ SX3 sx3_contr; TCutG *cut; TCutG *cut1; std::map, std::vector>> 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 verboseFit = true; // If true: print fit summary and chi² -const bool drawCanvases = false; // If false: canvases won't be drawn at all +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 = true; // If false: canvases won't be drawn at all void GainMatchSX3Front::Begin(TTree * /*tree*/) { @@ -102,7 +103,7 @@ void GainMatchSX3Front::Begin(TTree * /*tree*/) else backGainValid[id][bk] = false; } - + SX3 sx3_contr; } @@ -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,31 +293,30 @@ 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 xVals, yVals, exVals, eyVals; // Build data with fixed error for (size_t i = 0; i < udE.size(); ++i) { - double x = uE[i]; // front energy + double x = uE[i]; // front energy double y = dE[i]; // back energy 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; diff --git a/HistPlotter.h b/HistPlotter.h new file mode 100644 index 0000000..8bba450 --- /dev/null +++ b/HistPlotter.h @@ -0,0 +1,414 @@ +#ifndef HISTPLOTTER_H +#define HISTPLOTTER_H +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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 oMap; //!< Maps std::string to all TH1, TH2 objects in the class + std::unordered_map cutsMap; //!< Maps std::string to TCutG objects held by the class + std::set folderList; //!< List of all folder names used to nest objects + std::unordered_map 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> onedimcache; + std::unordered_map, std::vector>> 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(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(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(it->second)->FillN(onedimcache[it->first].size(), //size + onedimcache[it->first].data(), //array + std::vector(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(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(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 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(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(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(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(temp1D))); + onedimcache.insert(std::make_pair(name, std::vector())); + 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(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(temp1D))); + onedimcache.insert(std::make_pair(name, std::vector())); + onedimcache[name].reserve(16384); + if(foldername!="") { + if(folderList.find(foldername)==folderList.end()) { + folderList.insert(foldername); + } + foldersForObjects.insert(std::make_pair(static_cast(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(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(temp2D))); + twodimcache.insert(std::make_pair(name, std::make_pair(std::vector(),std::vector()))); + 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(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(temp2D))); + twodimcache.insert(std::make_pair(name, std::make_pair(std::vector(),std::vector()))); + 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(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(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(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