From aeaef18b41614b56a22c4189a677a394c14d4073 Mon Sep 17 00:00:00 2001
From: Eoin Clerkin <e.clerkin@ygsi.de>
Date: Fri, 25 Nov 2022 11:25:38 +0100
Subject: [PATCH] Removes buggy matbudget macros
All mcbm material budget map generating macros were found to have a
halfing of radlen bug and a missing volume bug. To prevent these macros
from being further used, they are removed.
---
macro/mcbm/matbudget_ana_mcbm_mvd.C | 165 -----------------------
macro/mcbm/matbudget_ana_mcbm_sts.C | 160 -----------------------
macro/mcbm/matbudget_ana_trd.C | 191 ---------------------------
macro/mcbm/matbudget_mc_mcbm_mvd.C | 195 ----------------------------
macro/mcbm/matbudget_mc_mcbm_sts.C | 194 ---------------------------
5 files changed, 905 deletions(-)
delete mode 100644 macro/mcbm/matbudget_ana_mcbm_mvd.C
delete mode 100644 macro/mcbm/matbudget_ana_mcbm_sts.C
delete mode 100644 macro/mcbm/matbudget_ana_trd.C
delete mode 100644 macro/mcbm/matbudget_mc_mcbm_mvd.C
delete mode 100644 macro/mcbm/matbudget_mc_mcbm_sts.C
diff --git a/macro/mcbm/matbudget_ana_mcbm_mvd.C b/macro/mcbm/matbudget_ana_mcbm_mvd.C
deleted file mode 100644
index e2459efa51..0000000000
--- a/macro/mcbm/matbudget_ana_mcbm_mvd.C
+++ /dev/null
@@ -1,165 +0,0 @@
-/* Copyright (C) 2015 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: David Emschermann [committer] */
-
-#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_mvd.C++("mvd_v18a")'
-
-Int_t matbudget_ana_mcbm_mvd(const char* inGeo, Int_t nEvents = 10000000)
-{
-
- // Input file (MC)
- TString geoVersion(inGeo);
- TString inFile = "data/matbudget." + geoVersion + ".mc.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 = 16; // maximal radius for histograms (for both x and y)
- const int rMax = 8; // 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; // STS
- name += i; // MVD
- 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 = posIn.Z()/10 - 3 + 0.5; // STS // suppose equidistant stations at 30-100 cm
- // int iStationOut = posOut.Z()/10 - 3 + 0.5; // STS
- int iStation = posIn.Z() / 5 - 1 + 0.5; // MVD
- int iStationOut = posOut.Z() / 5 - 1 + 0.5; // MVD
- 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"); // STS
- hStaRadLen[iStation]->SetAxisRange(0, 0.5, "Z"); // MVD
- 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_ana_mcbm_sts.C b/macro/mcbm/matbudget_ana_mcbm_sts.C
deleted file mode 100644
index 6bd0f2362b..0000000000
--- a/macro/mcbm/matbudget_ana_mcbm_sts.C
+++ /dev/null
@@ -1,160 +0,0 @@
-/* Copyright (C) 2015 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: David Emschermann [committer] */
-
-#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_sts(const char* inGeo, Int_t nEvents = 10000000)
-{
-
- // Input file (MC)
- TString geoVersion(inGeo);
- TString inFile = "data/matbudget." + geoVersion + ".mc.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 = 16; // 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 = posIn.Z() / 10 - 3 + 0.5; // suppose equidistant stations at 30-100 cm
- int iStationOut = posOut.Z() / 10 - 3 + 0.5;
- 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_ana_trd.C b/macro/mcbm/matbudget_ana_trd.C
deleted file mode 100644
index 89bc7de679..0000000000
--- a/macro/mcbm/matbudget_ana_trd.C
+++ /dev/null
@@ -1,191 +0,0 @@
-/* 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_trd(const TString geoVersion = "", const TString inFile = "matbudget", Int_t nEvents = -1)
-{
-
- // TRD stations Z and thickness
-
- const std::vector<double> trdZ = {441.6, 486.6, 531.6, 576.6};
-
- double trdZthick = 20; // station material is collected at trdZ +/- trdZthick
-
- const int nBins = 300; // number of bins in histograms (in both x and y)
- const int maxXY = 500; // +/- X and Y size for histograms (must be the same for X and Y)
-
- const int nStations = trdZ.size(); // number of stations
-
- for (int i = 1; i < nStations; i++) {
- if (trdZ[i] - trdZthick <= trdZ[i - 1] + trdZthick) {
- cerr << "Error:: TRD stations are not ordered in Z !!! " << endl;
- exit(-1);
- }
- }
-
- // Input file (MC)
- TFile* input = new TFile(inFile + ".tra.root");
- 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(-1);
- }
-
- // Output file (material maps)
- TString outFile = "trd_matbudget_" + geoVersion + ".root";
-
- // Input tree and branch
- TTree* tree = (TTree*) input->Get("cbmsim");
- TClonesArray* points = new TClonesArray("FairRadLenPoint");
- tree->SetBranchAddress("RadLen", &points);
-
- // Create output histograms
-
- TProfile2D* hStaRadLen[nStations];
- for (int i = 0; i < nStations; ++i) {
- TString name = "Radiation Thickness [%], Station";
- name += i + 1;
- hStaRadLen[i] = new TProfile2D(name, name, nBins, -maxXY, maxXY, nBins, -maxXY, maxXY);
- }
-
- // Event loop
- if (nEvents < 0 || nEvents > tree->GetEntriesFast()) { nEvents = tree->GetEntriesFast(); }
-
- vector<double> RadThick(nStations, 0.);
-
- for (Int_t event = 0; event < nEvents; event++) {
- if (0 == event % 100000) { std::cout << "*** Processing event " << event << std::endl; }
- tree->GetEntry(event);
-
- RadThick.assign(nStations, 0.);
-
- // For this implementation to be correct, there should be only one MCTrack per event.
-
- // Point loop. all points must xave the same (x,y) coordinates
-
- double x = 2 * maxXY;
- double y = 2 * maxXY;
- int trackID = -1;
-
- 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 = 0.5 * (posOut + posIn);
- const double xp = middle.X();
- const double yp = middle.Y();
- const double zp = middle.Z();
- const int idp = point->GetTrackID();
-
- if (iPoint == 0) {
- trackID = idp;
- x = xp;
- y = yp;
- }
-
- if (fabs(xp - x) > 1.e8 || fabs(yp - y) > 1.e8) {
- cerr << "Error: the track is not horisontal xy old: " << x << " " << y << " xy new: " << x << " " << y << endl;
- exit(-1);
- }
- if (idp != trackID) {
- cerr << "Error: more than one MC track in event" << endl;
- exit(-1);
- }
-
- // Material budget per unit length
- const double radThick = posDif.Mag() / point->GetRadLength();
-
- // Determine station number
- int iStationIn = -1;
- int iStationOut = -1;
-
- double zIn = posIn.Z();
- double zOut = posOut.Z();
-
- for (int i = 0; i < nStations; i++) {
- if (zIn > trdZ[i] - trdZthick) iStationIn = i;
- if (zOut > trdZ[i] - trdZthick) iStationOut = i;
- }
- if (zIn > trdZ[nStations - 1] + trdZthick) iStationIn = -1;
- if (zOut > trdZ[nStations - 1] + trdZthick) iStationOut = -1;
-
- if (iStationOut != iStationIn || iStationIn < 0) continue;
- assert(iStationIn >= nStations);
-
- RadThick[iStationIn] += radThick;
- }
-
- // Fill material budget map for each station in [%] of the radiation length
- for (int i = 0; i < nStations; ++i) {
- hStaRadLen[i]->Fill(x, y, RadThick[i] * 100);
- }
- } // 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, 5, "Z");
- hStaRadLen[iStation]->DrawCopy("colz");
- output->cd();
- 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_mc_mcbm_mvd.C b/macro/mcbm/matbudget_mc_mcbm_mvd.C
deleted file mode 100644
index b76c30ad83..0000000000
--- a/macro/mcbm/matbudget_mc_mcbm_mvd.C
+++ /dev/null
@@ -1,195 +0,0 @@
-/* Copyright (C) 2015 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: David Emschermann [committer] */
-
-// --------------------------------------------------------------------------
-//
-// Macro for transport simulation with the STS
-// to determine the material budget of the STS 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
-//
-// T. Balog, 11.05.2015
-//
-// --------------------------------------------------------------------------
-
-//void matbudget_mc_mcbm_sts(const char* inGeo, Int_t nEvents = 10000000) // +-50 x +-50
-void matbudget_mc_mcbm_mvd(const char* inGeo,
- Int_t nEvents = 1000000) // +-15 x +-15 = 0.09 ratio
-{
-
- // ========================================================================
- // Adjust this part according to your requirements
-
-
- // ----- Paths and file names --------------------------------------------
- TString geoVersion(inGeo);
- TString inDir = gSystem->Getenv("VMCWORKDIR");
- TString inFile = inDir + "/input/urqmd.ftn14";
- TString outDir = "data";
- TString outFile = outDir + "/matbudget." + geoVersion + ".mc.root";
- TString parFile = outDir + "/matbudget." + geoVersion + ".params.root";
-
- // ----- Geometries -----------------------------------------------------
- TString caveGeom = "cave.geo";
- TString targetGeom = "";
- TString pipeGeom = "";
- TString magnetGeom = "";
- TString mvdGeom = "mcbm/" + geoVersion + ".geo.root"; // "mcbm/mvd_v18a.geo.root";
- TString stsGeom = ""; // "mcbm/sts_" + geoVersion + ".geo.root";
- TString richGeom = "";
- TString trdGeom = "";
- TString tofGeom = "";
-
- // In general, the following parts need not be touched
- // ========================================================================
-
-
- // ---- Debug option -------------------------------------------------
- gDebug = 0;
- // ------------------------------------------------------------------------
-
-
- // ----- Timer --------------------------------------------------------
- TStopwatch timer;
- timer.Start();
- // ------------------------------------------------------------------------
-
-
- // ----- Create simulation run ----------------------------------------
- FairRunSim* run = new FairRunSim();
- run->SetName("TGeant3"); // Transport engine
- run->SetOutputFile(outFile); // Output file
- FairRuntimeDb* rtdb = run->GetRuntimeDb();
- // ------------------------------------------------------------------------
-
-
- // ----- Create media -------------------------------------------------
- run->SetMaterials("media.geo"); // Materials
- // ------------------------------------------------------------------------
-
-
- // ----- Create detectors and passive volumes -------------------------
- if (caveGeom != "") {
- FairModule* cave = new CbmCave("CAVE");
- cave->SetGeometryFileName(caveGeom);
- run->AddModule(cave);
- }
-
- if (pipeGeom != "") {
- FairModule* pipe = new CbmPipe("PIPE");
- pipe->SetGeometryFileName(pipeGeom);
- run->AddModule(pipe);
- }
-
- if (targetGeom != "") {
- FairModule* target = new CbmTarget("Target");
- target->SetGeometryFileName(targetGeom);
- run->AddModule(target);
- }
-
- if (magnetGeom != "") {
- FairModule* magnet = new CbmMagnet("MAGNET");
- magnet->SetGeometryFileName(magnetGeom);
- run->AddModule(magnet);
- }
-
- if (mvdGeom != "") {
- FairDetector* mvd = new CbmMvd("MVD", kTRUE);
- mvd->SetGeometryFileName(mvdGeom);
- run->AddModule(mvd);
- }
-
- if (stsGeom != "") {
- FairDetector* sts = new CbmStsMC(kTRUE);
- sts->SetGeometryFileName(stsGeom);
- run->AddModule(sts);
- }
-
- if (richGeom != "") {
- FairDetector* rich = new CbmRich("RICH", kTRUE);
- rich->SetGeometryFileName(richGeom);
- run->AddModule(rich);
- }
-
- if (trdGeom != "") {
- FairDetector* trd = new CbmTrd("TRD", kTRUE);
- trd->SetGeometryFileName(trdGeom);
- run->AddModule(trd);
- }
-
- if (tofGeom != "") {
- FairDetector* tof = new CbmTof("TOF", kTRUE);
- tof->SetGeometryFileName(tofGeom);
- run->AddModule(tof);
- }
-
- // ------------------------------------------------------------------------
-
-
- // ----- Create magnetic field ----------------------------------------
- // Zero field
- CbmFieldConst* magField = new CbmFieldConst();
- magField->SetField(0, 0, 0); // values are in kG
- magField->SetFieldRegion(-74, -39, -22, 74, 39, 160);
- run->SetField(magField);
-
- // ----- 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(-50.,-50.,50.,50.,0.); // STS SIS100
- // boxGen->SetBoxXYZ(-15.,-15.,15.,15.,0.); // STS mCBM
- boxGen->SetBoxXYZ(-7.5, -7.5, 7.5, 7.5, 0.); // MVD
- boxGen->SetPRange(0.1, 0.5);
- boxGen->SetThetaRange(0., 0.);
- boxGen->SetPhiRange(0., 360.);
-
- primGen->AddGenerator(boxGen);
-
- run->SetRadLenRegister(kTRUE);
- // ------------------------------------------------------------------------
-
- // ----- Run initialisation -------------------------------------------
- run->Init();
- // ------------------------------------------------------------------------
-
- // ----- Runtime database ---------------------------------------------
- CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar");
- fieldPar->SetParameters(magField);
- fieldPar->setChanged();
- fieldPar->setInputVersion(run->GetRunId(), 1);
- Bool_t kParameterMerged = kTRUE;
- FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
- parOut->open(parFile.Data());
- rtdb->setOutput(parOut);
- rtdb->saveOutput();
- rtdb->print();
- // ------------------------------------------------------------------------
-
-
- // ----- Start run ----------------------------------------------------
- run->Run(nEvents);
- // ------------------------------------------------------------------------
-
- // ----- Finish -------------------------------------------------------
- timer.Stop();
- Double_t rtime = timer.RealTime();
- Double_t ctime = timer.CpuTime();
- cout << endl << endl;
- cout << "Macro finished succesfully." << endl;
- cout << "Output file is " << outFile << endl;
- cout << "Parameter file is " << parFile << endl;
- cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << endl << endl;
- // ------------------------------------------------------------------------
-
- cout << " Test passed" << endl;
- cout << " All ok " << endl;
-}
diff --git a/macro/mcbm/matbudget_mc_mcbm_sts.C b/macro/mcbm/matbudget_mc_mcbm_sts.C
deleted file mode 100644
index 1fd9557a91..0000000000
--- a/macro/mcbm/matbudget_mc_mcbm_sts.C
+++ /dev/null
@@ -1,194 +0,0 @@
-/* Copyright (C) 2015 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: David Emschermann [committer] */
-
-// --------------------------------------------------------------------------
-//
-// Macro for transport simulation with the STS
-// to determine the material budget of the STS 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
-//
-// T. Balog, 11.05.2015
-//
-// --------------------------------------------------------------------------
-
-//void matbudget_mc_mcbm_sts(const char* inGeo, Int_t nEvents = 10000000) // +-50 x +-50
-void matbudget_mc_mcbm_sts(const char* inGeo,
- Int_t nEvents = 1000000) // +-15 x +-15 = 0.09 ratio
-{
-
- // ========================================================================
- // Adjust this part according to your requirements
-
-
- // ----- Paths and file names --------------------------------------------
- TString geoVersion(inGeo);
- TString inDir = gSystem->Getenv("VMCWORKDIR");
- TString inFile = inDir + "/input/urqmd.ftn14";
- TString outDir = "data";
- TString outFile = outDir + "/matbudget." + geoVersion + ".mc.root";
- TString parFile = outDir + "/matbudget." + geoVersion + ".params.root";
-
- // ----- Geometries -----------------------------------------------------
- TString caveGeom = "cave.geo";
- TString targetGeom = "";
- TString pipeGeom = "";
- TString magnetGeom = "";
- TString mvdGeom = "";
- TString stsGeom = "mcbm/" + geoVersion + ".geo.root";
- TString richGeom = "";
- TString trdGeom = "";
- TString tofGeom = "";
-
- // In general, the following parts need not be touched
- // ========================================================================
-
-
- // ---- Debug option -------------------------------------------------
- gDebug = 0;
- // ------------------------------------------------------------------------
-
-
- // ----- Timer --------------------------------------------------------
- TStopwatch timer;
- timer.Start();
- // ------------------------------------------------------------------------
-
-
- // ----- Create simulation run ----------------------------------------
- FairRunSim* run = new FairRunSim();
- run->SetName("TGeant3"); // Transport engine
- run->SetOutputFile(outFile); // Output file
- FairRuntimeDb* rtdb = run->GetRuntimeDb();
- // ------------------------------------------------------------------------
-
-
- // ----- Create media -------------------------------------------------
- run->SetMaterials("media.geo"); // Materials
- // ------------------------------------------------------------------------
-
-
- // ----- Create detectors and passive volumes -------------------------
- if (caveGeom != "") {
- FairModule* cave = new CbmCave("CAVE");
- cave->SetGeometryFileName(caveGeom);
- run->AddModule(cave);
- }
-
- if (pipeGeom != "") {
- FairModule* pipe = new CbmPipe("PIPE");
- pipe->SetGeometryFileName(pipeGeom);
- run->AddModule(pipe);
- }
-
- if (targetGeom != "") {
- FairModule* target = new CbmTarget("Target");
- target->SetGeometryFileName(targetGeom);
- run->AddModule(target);
- }
-
- if (magnetGeom != "") {
- FairModule* magnet = new CbmMagnet("MAGNET");
- magnet->SetGeometryFileName(magnetGeom);
- run->AddModule(magnet);
- }
-
- if (mvdGeom != "") {
- FairDetector* mvd = new CbmMvd("MVD", kTRUE);
- mvd->SetGeometryFileName(mvdGeom);
- run->AddModule(mvd);
- }
-
- if (stsGeom != "") {
- FairDetector* sts = new CbmStsMC(kTRUE);
- sts->SetGeometryFileName(stsGeom);
- run->AddModule(sts);
- }
-
- if (richGeom != "") {
- FairDetector* rich = new CbmRich("RICH", kTRUE);
- rich->SetGeometryFileName(richGeom);
- run->AddModule(rich);
- }
-
- if (trdGeom != "") {
- FairDetector* trd = new CbmTrd("TRD", kTRUE);
- trd->SetGeometryFileName(trdGeom);
- run->AddModule(trd);
- }
-
- if (tofGeom != "") {
- FairDetector* tof = new CbmTof("TOF", kTRUE);
- tof->SetGeometryFileName(tofGeom);
- run->AddModule(tof);
- }
-
- // ------------------------------------------------------------------------
-
-
- // ----- Create magnetic field ----------------------------------------
- // Zero field
- CbmFieldConst* magField = new CbmFieldConst();
- magField->SetField(0, 0, 0); // values are in kG
- magField->SetFieldRegion(-74, -39, -22, 74, 39, 160);
- run->SetField(magField);
-
- // ----- 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(-50.,-50.,50.,50.,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->SetRadLenRegister(kTRUE);
- // ------------------------------------------------------------------------
-
- // ----- Run initialisation -------------------------------------------
- run->Init();
- // ------------------------------------------------------------------------
-
- // ----- Runtime database ---------------------------------------------
- CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar");
- fieldPar->SetParameters(magField);
- fieldPar->setChanged();
- fieldPar->setInputVersion(run->GetRunId(), 1);
- Bool_t kParameterMerged = kTRUE;
- FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
- parOut->open(parFile.Data());
- rtdb->setOutput(parOut);
- rtdb->saveOutput();
- rtdb->print();
- // ------------------------------------------------------------------------
-
-
- // ----- Start run ----------------------------------------------------
- run->Run(nEvents);
- // ------------------------------------------------------------------------
-
- // ----- Finish -------------------------------------------------------
- timer.Stop();
- Double_t rtime = timer.RealTime();
- Double_t ctime = timer.CpuTime();
- cout << endl << endl;
- cout << "Macro finished succesfully." << endl;
- cout << "Output file is " << outFile << endl;
- cout << "Parameter file is " << parFile << endl;
- cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << endl << endl;
- // ------------------------------------------------------------------------
-
- cout << " Test passed" << endl;
- cout << " All ok " << endl;
-}
--
GitLab