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.

Sign up for GitHub

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

Jump to bottom

refactor: move material composition code to the left #1804

Merged
merged 3 commits into from
Jan 27, 2023
Merged

refactor: move material composition code to the left #1804

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5b6562a
move material composition code to the left
andiwand Jan 27, 2023
415faf2
delete after new
andiwand Jan 27, 2023
acf63c7
Merge branch 'main' into refactor/move-material-composition-code-to-t…
kodiakhq[bot] Jan 27, 2023
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
239 changes: 132 additions & 107 deletions Examples/Scripts/MaterialMapping/materialComposition.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,13 @@ struct MaterialHistograms {

MaterialHistograms() = default;

~MaterialHistograms() {
delete x0_vs_eta;
delete l0_vs_eta;
delete x0_vs_phi;
delete l0_vs_phi;
}

/// This fills the event into the histograms
/// and clears the cache accordingly
///
Expand Down Expand Up @@ -108,128 +115,146 @@ void materialComposition(const std::string& inFile, const std::string& treeName,
// Open the input file & get the tree
auto inputFile = TFile::Open(inFile.c_str());
auto inputTree = dynamic_cast<TTree*>(inputFile->Get(treeName.c_str()));
if (inputTree != nullptr) {
// Get the different atomic numbers
TCanvas* materialCanvas =
new TCanvas("materialCanvas", "Materials", 100, 100, 620, 400);
materialCanvas->cd();
// Draw all the atomic elements & get the histogram
inputTree->Draw("mat_A>>hA(100,0.5,100.5)");
TH1F* histA = dynamic_cast<TH1F*>(gDirectory->Get("hA"));
histA->Draw();

auto outputFile = TFile::Open(outFile.c_str(), "recreate");

float v_eta = 0;
float v_phi = 0;
std::vector<float>* stepLength = new std::vector<float>;
std::vector<float>* stepX0 = new std::vector<float>;
std::vector<float>* stepL0 = new std::vector<float>;
std::vector<float>* stepA = new std::vector<float>;

std::vector<float>* startX = new std::vector<float>;
std::vector<float>* startY = new std::vector<float>;
std::vector<float>* startZ = new std::vector<float>;

std::vector<float>* endX = new std::vector<float>;
std::vector<float>* endY = new std::vector<float>;
std::vector<float>* endZ = new std::vector<float>;

inputTree->SetBranchAddress("v_eta", &v_eta);
inputTree->SetBranchAddress("v_phi", &v_phi);
inputTree->SetBranchAddress("mat_step_length", &stepLength);

inputTree->SetBranchAddress("mat_X0", &stepX0);
inputTree->SetBranchAddress("mat_L0", &stepL0);
inputTree->SetBranchAddress("mat_A", &stepA);

inputTree->SetBranchAddress("mat_sx", &startX);
inputTree->SetBranchAddress("mat_sy", &startY);
inputTree->SetBranchAddress("mat_sz", &startZ);

inputTree->SetBranchAddress("mat_ex", &endX);
inputTree->SetBranchAddress("mat_ey", &endY);
inputTree->SetBranchAddress("mat_ez", &endZ);

// Loop over all entries ---------------
unsigned int entries = inputTree->GetEntries();
if (inputTree == nullptr) {
return;
}

// Get the different atomic numbers
TCanvas* materialCanvas =
new TCanvas("materialCanvas", "Materials", 100, 100, 620, 400);
materialCanvas->cd();
// Draw all the atomic elements & get the histogram
inputTree->Draw("mat_A>>hA(100,0.5,100.5)");
TH1F* histA = dynamic_cast<TH1F*>(gDirectory->Get("hA"));
histA->Draw();

auto outputFile = TFile::Open(outFile.c_str(), "recreate");

float v_eta = 0;
float v_phi = 0;
std::vector<float>* stepLength = new std::vector<float>;
andiwand marked this conversation as resolved.
Show resolved Hide resolved
std::vector<float>* stepX0 = new std::vector<float>;
std::vector<float>* stepL0 = new std::vector<float>;
std::vector<float>* stepA = new std::vector<float>;

std::vector<float>* startX = new std::vector<float>;
std::vector<float>* startY = new std::vector<float>;
std::vector<float>* startZ = new std::vector<float>;

std::vector<float>* endX = new std::vector<float>;
std::vector<float>* endY = new std::vector<float>;
std::vector<float>* endZ = new std::vector<float>;

inputTree->SetBranchAddress("v_eta", &v_eta);
inputTree->SetBranchAddress("v_phi", &v_phi);
inputTree->SetBranchAddress("mat_step_length", &stepLength);

inputTree->SetBranchAddress("mat_X0", &stepX0);
inputTree->SetBranchAddress("mat_L0", &stepL0);
inputTree->SetBranchAddress("mat_A", &stepA);

inputTree->SetBranchAddress("mat_sx", &startX);
inputTree->SetBranchAddress("mat_sy", &startY);
inputTree->SetBranchAddress("mat_sz", &startZ);

inputTree->SetBranchAddress("mat_ex", &endX);
inputTree->SetBranchAddress("mat_ey", &endY);
inputTree->SetBranchAddress("mat_ez", &endZ);

// Loop over all entries ---------------
unsigned int entries = inputTree->GetEntries();

#ifdef BOOST_AVAILABLE
std::cout << "*** Event Loop: " << std::endl;
progress_display event_loop_progress(entries * regions.size());
std::cout << "*** Event Loop: " << std::endl;
progress_display event_loop_progress(entries * regions.size());
#endif

// Loop of the regions
for (auto& r : regions) {
std::string rName = std::get<0>(r);
// Loop of the regions
for (auto& r : regions) {
std::string rName = std::get<0>(r);

// The material histograms ordered by atomic mass
std::map<unsigned int, MaterialHistograms> mCache;
// The material histograms ordered by atomic mass
std::map<unsigned int, MaterialHistograms> mCache;

// The material histograms for all
mCache[0] = MaterialHistograms(rName, 0, bins, eta);
for (unsigned int ib = 1; ib <= 100; ++ib) {
if (histA->GetBinContent(ib) > 0.) {
mCache[ib] = MaterialHistograms(rName, ib, bins, eta);
}
// The material histograms for all
mCache[0] = MaterialHistograms(rName, 0, bins, eta);
for (unsigned int ib = 1; ib <= 100; ++ib) {
if (histA->GetBinContent(ib) > 0.) {
mCache[ib] = MaterialHistograms(rName, ib, bins, eta);
}
}

for (unsigned int ie = 0; ie < entries; ++ie) {
// Get the entry
inputTree->GetEntry(ie);
for (unsigned int ie = 0; ie < entries; ++ie) {
// Get the entry
inputTree->GetEntry(ie);

#ifdef BOOST_AVAILABLE
++event_loop_progress;
++event_loop_progress;
#endif

// Accumulate the material per track
size_t steps = stepLength->size();
for (unsigned int is = 0; is < steps; ++is) {
float sX = startX->at(is);
float sY = startY->at(is);
float sZ = startZ->at(is);
float sR = sqrt(sX * sX + sY * sY);

float eX = endX->at(is);
float eY = endY->at(is);
float eZ = endZ->at(is);
float eR = sqrt(eX * eX + eY * eY);

float minR = std::get<1>(r);
float maxR = std::get<2>(r);
float minZ = std::get<3>(r);
float maxZ = std::get<4>(r);

if (minR > sR or minZ > sZ or maxR < eR or maxZ < eZ) {
continue;
}

float step = stepLength->at(is);
float X0 = stepX0->at(is);
float L0 = stepL0->at(is);

// The integral one
auto& all = mCache[0];
all.s_x0 += step / X0;
all.s_l0 += step / L0;

unsigned int sA = histA->FindBin(stepA->at(is));
// The current one
auto& current = mCache[sA];
current.s_x0 += step / X0;
current.s_l0 += step / L0;
}
// Fill the profiles and clear the cache
for (auto& [key, cache] : mCache) {
cache.fillAndClear(v_eta, v_phi);
// Accumulate the material per track
size_t steps = stepLength->size();
for (unsigned int is = 0; is < steps; ++is) {
float sX = startX->at(is);
float sY = startY->at(is);
float sZ = startZ->at(is);
float sR = sqrt(sX * sX + sY * sY);

float eX = endX->at(is);
float eY = endY->at(is);
float eZ = endZ->at(is);
float eR = sqrt(eX * eX + eY * eY);

float minR = std::get<1>(r);
float maxR = std::get<2>(r);
float minZ = std::get<3>(r);
float maxZ = std::get<4>(r);

if (minR > sR or minZ > sZ or maxR < eR or maxZ < eZ) {
continue;
}
}
// -----------------------------------

for (auto [key, cache] : mCache) {
cache.write();
float step = stepLength->at(is);
float X0 = stepX0->at(is);
float L0 = stepL0->at(is);

// The integral one
auto& all = mCache[0];
all.s_x0 += step / X0;
all.s_l0 += step / L0;

unsigned int sA = histA->FindBin(stepA->at(is));
// The current one
auto& current = mCache[sA];
current.s_x0 += step / X0;
current.s_l0 += step / L0;
}
// Fill the profiles and clear the cache
for (auto& [key, cache] : mCache) {
cache.fillAndClear(v_eta, v_phi);
}
}
outputFile->Close();
// -----------------------------------

for (auto [key, cache] : mCache) {
cache.write();
}
}

outputFile->Close();

delete stepLength;
delete stepX0;
delete stepL0;
delete stepA;

delete startX;
delete startY;
delete startZ;

delete endX;
delete endY;
delete endZ;

delete materialCanvas;
}