mCBM: add macro for TRD material map

macro/mcbm/matbudget_ana_mcbm_trd.C

+/* 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;
+}

macro/mcbm/matbudget_sim.C

+/* 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
+}