diff --git a/macro/mcbm/matbudget_ana_mcbm_trd.C b/macro/mcbm/matbudget_ana_mcbm_trd.C
new file mode 100644
index 0000000000000000000000000000000000000000..6bb26992a6cbc95a397e063348ef8775afea2393
--- /dev/null
+++ b/macro/mcbm/matbudget_ana_mcbm_trd.C
@@ -0,0 +1,181 @@
+/* Copyright (C) 2018-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Valentina Akishina [committer], David Emschermann */
+// --------------------------------------------------------------------------
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+//#include "FairRadLenPoint.h"
+
+#include "TCanvas.h"
+#include "TClonesArray.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "TProfile2D.h"
+#include "TROOT.h"
+#include "TString.h"
+#include "TStyle.h"
+#include "TSystem.h"
+#include "TTree.h"
+#include "TVector3.h"
+
+#include <iostream>
+#include <vector>
+using std::cout;
+using std::endl;
+using std::vector;
+#endif
+
+// need to set this root include path:
+// include path: -I/opt/cbm/fairsoft_jul15p1/installation/include/root -I/opt/cbm/fairroot_v-15.07-fairsoft_jul15p1/include -I/opt/cbm/fairsoft_jul15p1/installation/include
+// .include $SIMPATH/include
+// .include $FAIRROOTPATH/include
+// .x 'matbudget_ana_mcbm_sts.C++("sts_v18e")'
+
+Int_t matbudget_ana_mcbm_trd(const char* inGeo = " ", Int_t nEvents = 1000000)
+{
+
+ // Input file (MC)
+ TString geoVersion(inGeo);
+ TString inFile = "matbudget.tra.root";
+ TFile* input = new TFile(inFile);
+ if (!input) { // this is not checked in compiled mode
+ cout << "*** matbudget_ana: Input file " << inFile << " not found!\n"
+ << "Be sure to run matbudget_mc.C before for the respective " + geoVersion + " geometry!" << endl;
+ exit;
+ }
+
+ // Output file (material maps)
+ TString outFile = geoVersion + "_matbudget.root";
+
+ // Input tree and branch
+ TTree* tree = (TTree*) input->Get("cbmsim");
+ TClonesArray* points = new TClonesArray("FairRadLenPoint");
+ tree->SetBranchAddress("RadLen", &points);
+
+ // Create output histograms
+ TH1D* hisRadLen = new TH1D("hisRadLen", "Radiation Length", 1000, 0, 100);
+ // const int nStations = 8; // number of STS stations
+ // const int nBins = 1000; // number of bins in histograms (in both x and y)
+ // const int rMax = 55; // maximal radius for histograms (for both x and y)
+ const int nStations = 2; // number of STS stations
+ const int nBins = 300; // number of bins in histograms (in both x and y)
+ const int rMax = 100; // maximal radius for histograms (for both x and y)
+ TProfile2D* hStaRadLen[nStations];
+ for (int i = 0; i < nStations; ++i) {
+ TString name = "Radiation Thickness [%],";
+ name += " Station";
+ name += i + 1;
+ hStaRadLen[i] = new TProfile2D(name, name, nBins, -rMax, rMax, nBins, -rMax, rMax);
+ }
+
+ // Auxiliary variables
+ TVector3 vecs, veco;
+ std::map<int, int> trackHitMap;
+
+
+ // Event loop
+ int firstEvent = 0;
+ for (Int_t event = firstEvent; event < (firstEvent + nEvents) && event < tree->GetEntriesFast(); event++) {
+ tree->GetEntry(event);
+ if (0 == event % 100000) cout << "*** Processing event " << event << endl;
+
+ const int nTracks = 1;
+ std::vector<double> RadLengthOnTrack(nTracks, 0.0); //trackID, vector with points on track
+ std::vector<double> TrackLength(nTracks, 0.0); //trackID, vector with points on track
+
+ vector<double> RadThick(nStations, 0);
+ double x, y;
+
+ // For this implementation to be correct, there should be only one MCTrack per event.
+
+ // Point loop
+ for (Int_t iPoint = 0; iPoint < points->GetEntriesFast(); iPoint++) {
+ FairRadLenPoint* point = (FairRadLenPoint*) points->At(iPoint);
+
+ // Get track position at entry and exit of material
+ TVector3 posIn, posOut, posDif;
+ posIn = point->GetPosition();
+ posOut = point->GetPositionOut();
+ posDif = 0.5 * (posOut - posIn);
+
+ // Midpoint between in and out
+ const TVector3 middle = (posOut + posIn);
+ x = middle.X() / 2;
+ y = middle.Y() / 2;
+ const double z = middle.Z() / 2;
+
+ // Material budget per unit length
+ const double radThick = posDif.Mag() / point->GetRadLength();
+ RadLengthOnTrack[point->GetTrackID()] += radThick;
+ TrackLength[point->GetTrackID()] += posDif.Mag();
+
+
+ // Determine station number
+ int iStation = -1; //posIn.Z() / 10 - 3 + 0.5; // suppose equidistant stations at 30-100 cm
+ int iStationOut = -1; //posOut.Z() / 10 - 3 + 0.5;
+
+
+ float zIn = posIn.Z();
+ int ista = -1;
+
+ if (zIn > 110) ista = 0;
+ if (zIn > 147) ista = 1;
+ if (zIn > 160) ista = -1;
+
+ iStation = ista;
+
+ zIn = posOut.Z();
+ ista = -1;
+
+ if (zIn > 110) ista = 0;
+ if (zIn > 147) ista = 1;
+ if (zIn > 160) ista = -1;
+ iStationOut = ista;
+
+ if (iStationOut != iStation) continue;
+ if (iStation >= nStations || iStation < 0) continue;
+ RadThick[iStation] += radThick;
+ }
+
+ // Fill material budget map for each station
+ for (int i = 0; i < nStations; ++i)
+ hStaRadLen[i]->Fill(x, y, RadThick[i] * 100);
+
+ for (int k = 0; k < RadLengthOnTrack.size(); k++)
+ if (RadLengthOnTrack[k] > 0) hisRadLen->Fill(RadLengthOnTrack[k]);
+
+ } // event loop
+
+
+ // Plotting the results
+ TCanvas* can1 = new TCanvas();
+ // can1->Divide(nStations/2,2);
+ can1->Divide(2, nStations / 2);
+ gStyle->SetPalette(1);
+ gStyle->SetOptStat(0);
+
+ // Open output file
+ TFile* output = new TFile(outFile, "RECREATE");
+ output->cd();
+
+ for (int iStation = 0; iStation < nStations; iStation++) {
+ can1->cd(iStation + 1);
+ hStaRadLen[iStation]->GetXaxis()->SetTitle("x [cm]");
+ hStaRadLen[iStation]->GetYaxis()->SetTitle("y [cm]");
+ //hStaRadLen[i]->GetZaxis()->SetTitle("radiation thickness [%]");
+ hStaRadLen[iStation]->SetAxisRange(0, 2, "Z");
+ hStaRadLen[iStation]->Draw("colz");
+ hStaRadLen[iStation]->Write();
+ }
+
+ // Plot file
+ TString plotFile = geoVersion + "_matbudget.png";
+ can1->SaveAs(plotFile);
+
+ // Close files
+ input->Close();
+ output->Close();
+ cout << "Material budget maps written to " << outFile << endl;
+
+ return 0;
+}
diff --git a/macro/mcbm/matbudget_sim.C b/macro/mcbm/matbudget_sim.C
new file mode 100644
index 0000000000000000000000000000000000000000..1b5675909c6f6eed66b39dcc80cef6ea6e2b9aa4
--- /dev/null
+++ b/macro/mcbm/matbudget_sim.C
@@ -0,0 +1,162 @@
+/* Copyright (C) 2018-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Valentina Akishina [committer], Volker Friese, David Emschermann */
+
+// --------------------------------------------------------------------------
+//
+// V. Friese 15/07/2018
+//
+// The output file will be named [output].tra.root.
+// A parameter file [output].par.root will be created.
+// The geometry (TGeoManager) will be written into trd.geo.root.
+//
+// T. Balog, 11.05.2015
+//
+// Macro for transport simulation with the TRD
+// to determine the material budget of the TRD stations.
+// ROOTinos parallel to the z axis will be transported,
+// using the RadLenRegister functionality of FairRunSim.
+// From the output MC file, the material budget can be
+// determined with the macro matbudget_ana.C
+// --------------------------------------------------------------------------
+
+void SetTrack(CbmTransport*, Double_t, Int_t, Double_t, Double_t, Double_t);
+
+void matbudget_sim(const char* inGeo = "", const char* output = "matbudget", const char* inputFile = "",
+ const char* setupName = "mcbm_beam_2022_07", Int_t nEvents = 1000000)
+{
+ // --- Define the beam angle ----------------------------------------------
+ Double_t beamRotY = 25.;
+ // ------------------------------------------------------------------------
+
+ // --- Define the target geometry -----------------------------------------
+ //
+ // The target is not part of the setup, since one and the same setup can
+ // and will be used with different targets.
+ // The target is constructed as a tube in z direction with the specified
+ // diameter (in x and y) and thickness (in z). It will be placed at the
+ // specified position as daughter volume of the volume present there. It is
+ // in the responsibility of the user that no overlaps or extrusions are
+ // created by the placement of the target.
+ //
+ TString targetElement = "Gold";
+ Double_t targetPosX = 0.; // target x position in global c.s. [cm]
+ Double_t targetPosY = 0.; // target y position in global c.s. [cm]
+ Double_t targetPosZ = 0.; // target z position in global c.s. [cm]
+
+ // Double_t targetThickness = 0.1; // full thickness in cm
+ // Double_t targetDiameter = 0.5; // diameter in cm
+ // Double_t targetRotY = 25.; // target rotation angle around the y axis [deg]
+
+ Double_t targetThickness = 0.025; // mCBM thin gold target 0.25 mm = 0.025 cm thickness
+ Double_t targetDiameter = 1.5; // mCBM target width 15 mm = 1.5 cm
+ // Double_t targetDiameter = 0.1; // set small target for window acceptance plots
+ Double_t targetRotY = beamRotY; // target rotation angle around the y axis [deg]
+ // ------------------------------------------------------------------------
+
+ // --- Logger settings ----------------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("VERYHIGH");
+ // ------------------------------------------------------------------------
+
+ // ----- Environment --------------------------------------------------
+ TString myName = "mcbm_transport"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+
+ // ----- In- and output file names ------------------------------------
+ TString dataset(output);
+ TString outFile = dataset + ".tra.root";
+ TString parFile = dataset + ".par.root";
+ TString geoFile = dataset + ".geo.root";
+ std::cout << std::endl;
+}
+// ------------------------------------------------------------------------
+
+
+// ----- Timer --------------------------------------------------------
+TStopwatch timer;
+timer.Start();
+// ------------------------------------------------------------------------
+
+
+// --- Transport run ----------------------------------------------------
+CbmTransport run;
+
+
+// ----- Create PrimaryGenerator --------------------------------------
+FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
+// run.SetGenerator(primGen);
+
+// --- Primary particles
+// Generated are ROOTinos (PDG=0) in z direction (one per event),
+// starting at z = 0.
+// The starting points in x and y are chosen such as to illuminate the STS.
+FairBoxGenerator* boxGen = new FairBoxGenerator(0, 1);
+boxGen->SetBoxXYZ(-100., -100., 100., 100., 0.); // STS SIS100
+//boxGen->SetBoxXYZ(-15., -15., 15., 15., 0.); // STS mCBM
+boxGen->SetPRange(0.1, 0.5);
+boxGen->SetThetaRange(0., 0.);
+boxGen->SetPhiRange(0., 360.);
+
+primGen->AddGenerator(boxGen);
+
+run.RegisterRadLength(kTRUE);
+
+// comment the following line to remove target interaction
+run.AddInput(boxGen);
+run.SetOutFileName(outFile);
+run.SetParFileName(parFile);
+run.SetGeoFileName(geoFile);
+run.LoadSetup(setupName);
+run.SetField(new CbmFieldConst());
+// run.SetTarget(targetElement, targetThickness, targetDiameter, targetPosX, targetPosY, targetPosZ,
+// targetRotY * TMath::DegToRad());
+//run.SetBeamPosition(0., 0., 0.1, 0.1); // Beam width 1 mm is assumed
+// run.SetBeamAngle(beamRotY * TMath::DegToRad(), 0.);
+//run.StoreTrajectories();
+run.Run(nEvents);
+// ------------------------------------------------------------------------
+
+
+// ----- Finish -------------------------------------------------------
+timer.Stop();
+Double_t rtime = timer.RealTime();
+Double_t ctime = timer.CpuTime();
+std::cout << std::endl << std::endl;
+std::cout << "Macro finished successfully." << std::endl;
+std::cout << "Output file is " << outFile << std::endl;
+std::cout << "Parameter file is " << parFile << std::endl;
+std::cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << std::endl << std::endl;
+// ------------------------------------------------------------------------
+
+
+// ----- Resource monitoring ------------------------------------------
+FairSystemInfo sysInfo;
+Float_t maxMemory = sysInfo.GetMaxMemory();
+std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+std::cout << maxMemory;
+std::cout << "</DartMeasurement>" << std::endl;
+
+Float_t cpuUsage = ctime / rtime;
+std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+std::cout << cpuUsage;
+std::cout << "</DartMeasurement>" << std::endl;
+
+
+std::cout << " Test passed" << std::endl;
+std::cout << " All ok " << std::endl;
+// ------------------------------------------------------------------------
+}
+
+
+void SetTrack(CbmTransport* run, Double_t beamRotY, Int_t pdgid, Double_t x, Double_t y, Double_t z)
+{
+ TVector3 v;
+ v.SetXYZ(x, y, z);
+ v.RotateY(-beamRotY * acos(-1.) / 180.);
+ cout << "X " << v.X() << " Y " << v.Y() << " Z " << v.Z() << endl;
+
+ run->AddInput(new FairParticleGenerator(pdgid, 1, v.X(), v.Y(), v.Z())); // single electron along beam axis
+}