Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Response matrix macro #313

Merged
merged 6 commits into from
Oct 7, 2022
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 70 additions & 0 deletions macros/REST_GenerateResponseMatrix.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
//*******************************************************************************************************
//*** Description: This macro receives as input two variable names that must be present inside the
//*** analysis tree. It creates a TH2D histogram using those variables. The histogram limits and range
//*** may be specified in the last argument, following the ROOT `Draw` method, range definition format.
//***
//*** This method is typically used to generate the response matrix of detector simulations using
//*** Geant4, although it could be used to generate a matrix in binary format of a TH2 object defined
//*** by any two analysis tree variables.
//*** --------------
//*** The output file will be a binary file containing the table, this table could be read later on
//*** using the method TRestTools::ReadBinaryTable.
//*** --------------
//*** Usage: restManager GenerateResponseMatrix /full/path/file.root [varX] [varY] [range]
//***
//*** Input arguments:
//*** - `varX` and `varY` : two variables inside the analysis tree.
//*** - `range` : If given, it defines the histogram limits and the binning. It is defined as:
//*** (nBinsX, Xlow, Xhigh, nBinsY, yLow, yHigh)
//***
//**********************************************************************************************************
Int_t REST_GenerateResponseMatrix(std::string fname, std::string varX = "g4Ana_totalEdep",
std::string varY = "g4Ana_energyPrimary",
std::string range = "(150,0,15,150,0,15)") {
TRestRun run(fname);

TRestAnalysisTree* aTree = run.GetAnalysisTree();

TRestGeant4Metadata* g4Md = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");

std::string drawCommand = varY + ":" + varX;
if (range != "") drawCommand += ">>response" + range;

aTree->Draw((TString)drawCommand);

TH2D* h = (TH2D*)aTree->GetHistogram();

/// We renormalize the values so that the values will be given
/// on the units of X and Y axis.
Double_t lowXValue = h->GetXaxis()->GetBinLowEdge(1);
Double_t highXValue = h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
Double_t normX = (highXValue - lowXValue) / h->GetNbinsX();

Double_t lowYValue = h->GetYaxis()->GetBinLowEdge(1);
Double_t highYValue = h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());
Double_t normY = (highYValue - lowYValue) / h->GetNbinsY();

std::vector<std::vector<Float_t> > responseData;
for (int n = 1; n <= h->GetNbinsX(); n++) {
std::vector<Float_t> primaryResponse;
for (int m = 1; m <= h->GetNbinsY(); m++) {
Double_t value = h->GetBinContent(n, m) / normX / normY / g4Md->GetNumberOfEvents();
primaryResponse.push_back(value);
}
responseData.push_back(primaryResponse);
}

std::string output_fname =
(string)run.GetRunTag() + ".N" + REST_StringHelper::IntegerToString(responseData[0].size()) + "f";

std::cout << "Writting output binary file: " << output_fname << std::endl;

TRestTools::ExportBinaryTable(output_fname, responseData);

Double_t efficiency = h->Integral() / g4Md->GetNumberOfEvents();

std::cout << "Overall efficiency : " << efficiency << std::endl;
std::cout << "Number of primaries: " << g4Md->GetNumberOfEvents() << std::endl;

return 0;
}
24 changes: 24 additions & 0 deletions macros/REST_ListMacros.C
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,34 @@ using namespace std;

#include <TString.h>

//*******************************************************************************************************
//***
//*** Lists all the official macros together with its documentation
//***
//*******************************************************************************************************
Int_t REST_ListMacros() {
string macrosPath = (string)getenv("REST_PATH") + "/macros";
vector<string> directories = TRestTools::GetSubdirectories(macrosPath);

cout << "Directory : " << macrosPath << endl;
vector<string> main_files = TRestTools::GetFilesMatchingPattern(macrosPath + "/REST_*.C");
for (int m = 0; m < main_files.size(); m++) {
cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
cout << " ++ "
<< " Macro : " << TRestTools::SeparatePathAndName(main_files[m]).second << endl;
std::ifstream t(main_files[m]);
std::string str((std::istreambuf_iterator<char>(t)), std::istreambuf_iterator<char>());

std::vector<string> lines = REST_StringHelper::Split(str, "\n");

cout << " ++ " << endl;
for (int l = 0; l < lines.size(); l++)
if (lines[l].find("//*** ") != string::npos)
cout << " ++ " << lines[l].substr(6, -1) << endl;
cout << " ++ " << endl;
}
cout << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;

for (int n = 0; n < directories.size(); n++) {
cout << "Directory : " << directories[n] << endl;
if (directories[n].find("pipeline") != string::npos) continue;
Expand Down
6 changes: 4 additions & 2 deletions source/framework/tools/src/TRestTools.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -668,7 +668,7 @@ bool TRestTools::fileExists(const string& filename) { return std::filesystem::ex
///////////////////////////////////////////////
/// \brief Returns true if the **filename** has *.root* extension.
///
bool TRestTools::isRootFile(const string& filename) { return GetFileNameExtension(filename) == ".root"; }
bool TRestTools::isRootFile(const string& filename) { return GetFileNameExtension(filename) == "root"; }

///////////////////////////////////////////////
/// \brief Returns true if **filename** is an *http* address.
Expand Down Expand Up @@ -729,7 +729,9 @@ std::pair<string, string> TRestTools::SeparatePathAndName(const string& fullname
/// Input: "/home/jgalan/abc.txt" Output: "txt"
///
string TRestTools::GetFileNameExtension(const string& fullname) {
return filesystem::path(fullname).extension().string();
string extension = filesystem::path(fullname).extension().string();
if (extension.size() > 1) return extension.substr(1);
return extension;
}

///////////////////////////////////////////////
Expand Down