tracking_detectors/scripts/test_matscan.cxx

0001 // material scan of a cylinder volume for athena
0002 // the scan starts at 0 for now
0003 // based on https://dd4hep.web.cern.ch/dd4hep/reference/MaterialScan_8cpp_source.html
0004 // dd4hep also provides command line utilities materialBudget and materialScan
0005 //
0006 // Shujie Li, Aug 2021
0007
0008 R__LOAD_LIBRARY(libDDCore.so)
0009 R__LOAD_LIBRARY(libDDG4.so)
0010 R__LOAD_LIBRARY(libDDG4IO.so)
0011 #include "DD4hep/Detector.h"
0012 #include "DD4hep/DetElement.h"
0013 #include "DD4hep/Objects.h"
0014 #include "DD4hep/Detector.h"
0015 #include "DDG4/Geant4Data.h"
0016 #include "DDRec/CellIDPositionConverter.h"
0017 #include "DDRec/SurfaceManager.h"
0018 #include "DDRec/Surface.h"
0019
0020 #include "TCanvas.h"
0021 #include "TChain.h"
0022 #include "TGeoMedium.h"
0023 #include "TGeoManager.h"
0024 #include "DDRec/MaterialScan.h"
0025 #include "DDRec/MaterialManager.h"
0026 #include "DD4hep/Detector.h"
0027 #include "DD4hep/Printout.h"
0028 #include "fmt/core.h"
0029
0030 #include <iostream>
0031 #include <fstream>
0032
0033 using namespace dd4hep;
0034 using namespace dd4hep::rec;
0035
0036   // for silicon tracker, the outer barrel rmax ~45, outer endcap zmax~121
0037 void test_matscan(const char* compact = "athena.xml",double zmax=130, double rmax = 61){
0038
0039   dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0040   detector.fromCompact(compact);
0041   MaterialScan matscan(detector);
0042   TString detname = "all";
0043   fmt::print("\n");
0044   fmt::print("All detector subsystem names:\n");
0045   for(const auto&  d : detector.detectors() ) {
0046     fmt::print("  {}\n", d.first);
0047   }
0048   // to do: material scan of given dtectors.
0049
0050   // TString det_list[14]={"TrackerBarrel_Inner","TrackerBarrel_Outer","TrackerEndcapN_Inner","TrackerEndcapN_Outer","TrackerEndcapP_Inner","TrackerEndcapP_Outer","TrackerSubAssembly_Inner","TrackerSubAssembly_Outer","VertexBarrel","VertexBarrelSubAssembly","VertexEndcapN","VertexEndcapP","VertexEndcapSubAssembly","cb_DIRC"};
0051   // for (int dd=0;dd<14;dd++){
0052   // TString detname = det_list[dd];
0053   // matscan.setDetector(detname.Data());
0054
0055   const char* fmt1 = "%8.3f %8.3f %8.3f %3d %-20s %3.0f  %8.4f %11.4f %11.4f %11.4f %11.4f %8.3f %8.3f %8.3f\n";
0056
0057   // the beam vacuum is missing if starting from origin and ending at negative z. need to check later
0058   double x0=0,y0=0,z0=0.001,x1,y1,z1;  // cm
0059   double epsilon=1e-4; // (mm) default 1e-4: Materials with a thickness smaller than epsilon (default 1e-4=1mu
0060
0061   const double DegToRad=3.141592653589793/180.0;
0062   double phi,dphi=0.5;    // Degree
0063   double dz=0.5, dr=0.5;    // cm
0064
0065   vector<double> lx,ly,lz;
0066   // prepare coord. for barrel
0067   for(z1=-zmax;z1<=zmax;z1+=dz){
0068     for(phi=0.;phi<=360;phi+=dphi){
0069       x1 = cos(phi*DegToRad)*rmax;
0070       y1 = sin(phi*DegToRad)*rmax;
0071       lx.push_back(x1);
0072       ly.push_back(y1);
0073       lz.push_back(z1);
0074     }
0075   }
0076   // prepare coord. for endcaps
0077   for (double r=1; r<=rmax; r+=dr){
0078     for(phi=0.;phi<=360;phi+=dphi){
0079       x1 = cos(phi*DegToRad)*r;
0080       y1 = sin(phi*DegToRad)*r;
0081       lx.push_back(x1);
0082       ly.push_back(y1);
0083       lz.push_back(zmax);
0084       lx.push_back(x1);
0085       ly.push_back(y1);
0086       lz.push_back(-zmax);
0087     }
0088   }
0089
0090   // loop over coord ponits for material scan
0091   long nn = lx.size();
0092   TString fname = Form("%s_z%g_r%g.dat",detname.Data(),zmax,rmax);
0093   FILE * pFile;
0094   pFile = fopen (fname,"w");
0095
0096   for(int ii=0;ii<nn;ii++){
0097     x1=lx[ii]; y1=ly[ii]; z1=lz[ii];
0098     if(ii%5000==0)
0099       cout<<ii<<"/"<<nn<<":"<<x1<<" "<<y1<<" "<<z1<<endl;
0100
0101     // if ((x1*x1+y1*y1)>(60*60))
0102       // continue;
0103     Vector3D p0(x0, y0, z0), p1(x1, y1, z1);
0104     Vector3D end, direction;
0105     direction = (p1-p0).unit();
0106     const auto& placements = matscan.scan(x0,y0,z0,x1,y1,z1,epsilon);
0107     // matscan.print(x0,y0,z0,x1,y1,z1,epsilon);
0108
0109     double sum_x0 = 0;
0110     double sum_lambda = 0;
0111     double path_length = 0, total_length = 0;
0112
0113
0114     for (unsigned i=0;i<placements.size();i++){
0115
0116       TGeoMaterial* mat=placements[i].first->GetMaterial();
0117       double length = placements[i].second;
0118       double nx0     = length / mat->GetRadLen();
0119       double nLambda = length / mat->GetIntLen();
0120       sum_x0        += nx0;
0121       sum_lambda    += nLambda;
0122       path_length   += length;
0123       total_length  += length;
0124       end = p0 + total_length * direction;
0125
0126
0127       fprintf(pFile, fmt1,x1, y1, z1, i+1, mat->GetName(), mat->GetZ(),
0128                 mat->GetDensity(), mat->GetRadLen(),
0129                 length, path_length, sum_x0,end[0], end[1], end[2]);
0130
0131
0132
0133     }
0134     // cout<<detname<<"  "<<x1<<","<<y1<<","<<z1<<endl;
0135     // cout<<x1<<","<<y1<<","<<z1<<": "<<placements.size()<<"  "<<sum_x0<<"  "<<total_length<<endl;
0136   }
0137   fclose (pFile);
0138 }
0139