diff --git a/Examples/Scripts/MaterialMapping/materialComposition.C b/Examples/Scripts/MaterialMapping/materialComposition.C index 7b77cda1166..beff638d09f 100644 --- a/Examples/Scripts/MaterialMapping/materialComposition.C +++ b/Examples/Scripts/MaterialMapping/materialComposition.C @@ -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 /// @@ -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(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(gDirectory->Get("hA")); - histA->Draw(); - - auto outputFile = TFile::Open(outFile.c_str(), "recreate"); - - float v_eta = 0; - float v_phi = 0; - std::vector* stepLength = new std::vector; - std::vector* stepX0 = new std::vector; - std::vector* stepL0 = new std::vector; - std::vector* stepA = new std::vector; - - std::vector* startX = new std::vector; - std::vector* startY = new std::vector; - std::vector* startZ = new std::vector; - - std::vector* endX = new std::vector; - std::vector* endY = new std::vector; - std::vector* endZ = new std::vector; - - 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(gDirectory->Get("hA")); + histA->Draw(); + + auto outputFile = TFile::Open(outFile.c_str(), "recreate"); + + float v_eta = 0; + float v_phi = 0; + std::vector* stepLength = new std::vector; + std::vector* stepX0 = new std::vector; + std::vector* stepL0 = new std::vector; + std::vector* stepA = new std::vector; + + std::vector* startX = new std::vector; + std::vector* startY = new std::vector; + std::vector* startZ = new std::vector; + + std::vector* endX = new std::vector; + std::vector* endY = new std::vector; + std::vector* endZ = new std::vector; + + 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 mCache; + // The material histograms ordered by atomic mass + std::map 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; }