UtilityApps/src/materialBudget.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 //==========================================================================
0011 //
0012 //  Compute the  material budget in terms of integrated radiation and ineraction
0013 //  lengths as seen from the IP for various subdetector regions given in
0014 //  a simple steering file. Creates a root file w/ histograms and dumps numbers
0015 //  the screen.
0016 //
0017 //  Author     : F.Gaede, DESY
0018 //
0019 //==========================================================================
0020
0021 // Framework include files
0022 #include <DD4hep/DetType.h>
0023 #include <DD4hep/Printout.h>
0024 #include <DD4hep/Detector.h>
0025 #include <DDRec/MaterialManager.h>
0026
0027 #include <TFile.h>
0028 #include <TH1F.h>
0029 #include <TError.h>
0030
0031 #include <cerrno>
0032 #include <fstream>
0033 #include <climits>
0034
0035 #include "main.h"
0036
0037 using namespace dd4hep;
0038 using namespace dd4hep::rec;
0039
0040 void dumpExampleSteering() ;
0041
0042 namespace {
0043
0044   struct SDetHelper{
0045     std::string name { };
0046     double r0 { 0.0 };
0047     double r1 { 0.0 };
0048     double z0 { 0.0 };
0049     double z1 { 0.0 };
0050     TH1F*  hx { nullptr };
0051     TH1F*  hl { nullptr };
0052   };
0053
0054   /// comput a point on a cylinder (r,z) for a given direction (theta,phi) from the IP
0055   Vector3D pointOnCylinder( double theta, double r, double z, double phi){
0056
0057     double  theta0 = atan2( r, z ) ;
0058
0059     Vector3D v = (   theta > theta0      ?
0060                      Vector3D(       r           , phi , r / tan( theta ) , Vector3D::cylindrical )  :
0061                      Vector3D( z * tan( theta )  , phi ,      z           , Vector3D::cylindrical )  ) ;
0062     return v ;
0063   }
0064
0065 }
0066
0067 int main_wrapper(int argc, char** argv)   {
0068
0069   struct Handler  {
0070     Handler() { SetErrorHandler(Handler::print); }
0071
0072     static void print(int level, Bool_t abort, const char *location, const char *msg)  {
0073       if ( level > kInfo || abort ) ::printf("%s: %s\n", location, msg);
0074     }
0075
0076     static void usage()  {
0077       std::cout << " usage: materialBudget compact.xml steering.txt" << std::endl
0078                 << "     -> create histograms with the material budget as seen from the IP within fixed ranges of (rmin, rmax, zmin. zmax)" << std::endl
0079                 << "        see example steering file for details ..." << std::endl
0080                 << "  -x   : dump example steering file " << std::endl
0081                 << std::endl;
0082       exit(1);
0083     }
0084
0085   } _handler;
0086
0087
0088   if( argc == 2 ){
0089     std::string dummy = argv[1] ;
0090
0091     if( dummy == "-x" )
0092       dumpExampleSteering() ;
0093     else
0094       Handler::usage();
0095   }
0096
0097   if( argc != 3 )
0098     Handler::usage();
0099
0100
0101   // ================== parse steering file ========================================================
0102   std::string compactFile  =  argv[1];
0103   std::string steeringFile =  argv[2];
0104   int nbins = 90 ;
0105   double phi0 = M_PI / 2. ;
0106   double thetaMin = 0. ;
0107   double thetaMax = 90. ;
0108   double etaMin = 0. ;
0109   double etaMax = -1. ;
0110   std::string outFileName("material_budget.root") ;
0111   std::vector<SDetHelper> subdets ;
0112
0113
0114   std::ifstream infile(steeringFile );
0115   std::string line;
0116
0117   while (std::getline(infile, line))  {
0118     // ignore empty lines and comments
0119     if( line.empty() || line.find("#") == 0 )
0120       continue ;
0121
0122     std::istringstream iss(line);
0123     std::string token ;
0124     iss >> token ;
0125
0126     if( token == "nbins" ){
0127       iss >> nbins ;
0128     }
0129     else if( token == "thetaMin" ){
0130       iss >> thetaMin ;
0131     }
0132     else if( token == "thetaMax" ){
0133       iss >> thetaMax ;
0134     }
0135     else if( token == "etaMin" ){
0136       iss >> etaMin ;
0137     }
0138     else if( token == "etaMax" ){
0139       iss >> etaMax ;
0140     }
0141     else if( token == "rootfile" ){
0142       iss >> outFileName ;
0143     }
0144     else if( token == "phi" ){
0145       iss >> phi0 ;
0146       phi0 = phi0 / 180. * M_PI ;
0147     }
0148     else if( token == "subdet" ){
0149       SDetHelper det;
0150       iss >>  det.name >> det.r0 >> det.z0 >> det.r1 >> det.z1 ;
0151       subdets.emplace_back( det );
0152     }
0153
0154     if ( iss.fail() ){
0155       std::cout << " ERROR parsing line : " << line << std::endl ;
0156       ::exit(EINVAL);
0157     }
0158   }
0159
0160   if ( nbins <= 0 )    {
0161     std::cout << "Invalid value for binning (nbins) " << nbins << std::endl;
0162     return EINVAL;
0163   }
0164
0165   // =================================================================================================
0166   setPrintLevel(WARNING);
0167
0168   Detector& description = Detector::getInstance();
0169   description.fromXML(compactFile ) ;
0170
0171   //----- open root file and book histograms
0172   TFile* rootFile = new TFile(outFileName.c_str(),"recreate");
0173   for( auto& det : subdets )  {
0174     std::string hxn(det.name), hxnn(det.name) ;
0175     std::string hln(det.name), hlnn(det.name) ;
0176
0177     if( etaMax > 0. )  {  // use eta
0178       hxn += "x0" ;
0179       hxnn += " integrated X0 vs eta" ;
0180       det.hx = new TH1F( hxn.c_str(), hxnn.c_str(), nbins, etaMin , etaMax ) ;
0181
0182       hln += "lambda" ;
0183       hlnn += " integrated int. lengths vs eta" ;
0184       det.hl = new TH1F( hln.c_str(), hlnn.c_str(), nbins, etaMin , etaMax ) ;
0185     }
0186     else {   // use polar angle
0187       hxn += "x0" ;
0188       hxnn += " integrated X0 vs -theta" ;
0189       det.hx = new TH1F( hxn.c_str(), hxnn.c_str(), nbins, -thetaMax , -thetaMin ) ;
0190
0191       hln += "lambda" ;
0192       hlnn += " integrated int. lengths vs -theta" ;
0193       det.hl = new TH1F( hln.c_str(), hlnn.c_str(), nbins, -thetaMax , -thetaMin ) ;
0194     }
0195   }
0196   //-------------------------
0197
0198
0199   Volume world  = description.world().volume() ;
0200
0201   MaterialManager matMgr( world ) ;
0202
0203   thetaMin = thetaMin / 180. * M_PI ;
0204   thetaMax = thetaMax / 180. * M_PI ;
0205   double dTheta = (thetaMax-thetaMin)/nbins; // bin size
0206   double dEta  = (etaMax-etaMin)/nbins ;
0207
0208   std::cout  << "====================================================================================================" << std::endl ;
0209   std::cout  << "theta:f/" ;
0210   for(auto& det : subdets){ std::cout  << det.name << "_x0:f/" << det.name << "_lam:f/" ; }
0211   std::cout  << std::endl ;
0212
0213   if ( nbins <= 0 || nbins > USHRT_MAX )  {
0214     std::cout << "Unreasonable number of bins: " << nbins << std::endl;
0215     ::exit(EINVAL);
0216   }
0217
0218   for(int i=0 ; i< nbins ;++i){
0219     double theta = ( etaMax > 0. ?  2. * atan ( exp ( - ( etaMin + (0.5+i)*dEta) ) ) : ( thetaMin + (0.5+i)*dTheta ) ) ;
0220     std::stringstream paramLine;
0221
0222     paramLine << std::scientific << theta << " " ;
0223     for( auto& det : subdets )  {
0224       Vector3D p0 = pointOnCylinder( theta, det.r0 , det.z0 , phi0  ) ;// double theta, double r, double z, double phi)
0225       Vector3D p1 = pointOnCylinder( theta, det.r1 , det.z1 , phi0  ) ;// double theta, double r, double z, double phi)
0226       const MaterialVec& materials = matMgr.materialsBetween(p0, p1);
0227       double sum_x0(0.), sum_lambda(0.);
0228       // double path_length(0.);
0229
0230       for( auto amat : materials )  {
0231         TGeoMaterial* mat =  amat.first->GetMaterial();
0232         double length = amat.second;
0233         double nx0 = length / mat->GetRadLen();
0234         sum_x0 += nx0;
0235         double nLambda = length / mat->GetIntLen();
0236         sum_lambda += nLambda;
0237         // path_length += length;
0238       }
0239
0240       double binX = ( etaMax > 0. ? (etaMin + (0.5+i)*dEta) : -theta/M_PI*180. ) ;
0241       det.hx->Fill( binX , sum_x0 ) ;
0242       det.hl->Fill( binX , sum_lambda ) ;
0243       paramLine  << std::scientific  << sum_x0 << "  " << sum_lambda << "  " ; // << path_length ;
0244     }
0245     std::cout << paramLine.str() << std::endl;
0246   }
0247   std::cout  << "====================================================================================================" << std::endl ;
0248   rootFile->Write();
0249   rootFile->Close();
0250   return 0;
0251 }
0252
0253 void dumpExampleSteering(){
0254   std::cout << "# Example steering file for materialBudget  (taken from ILD_l5_v02)" << std::endl ;
0255   std::cout <<  std::endl ;
0256   std::cout << "# root output file"  << std::endl ;
0257   std::cout << "rootfile material_budget.root" << std::endl ;
0258   std::cout << "# number of bins for polar angle (default 90)" << std::endl ;
0259   std::cout << "nbins 90" << std::endl ;
0260   std::cout <<  std::endl ;
0261   std::cout << "# use polar angle - specify minimum and maximum theta value" << std::endl ;
0262   std::cout << "# thetaMin 0." << std::endl ;
0263   std::cout << "# thetaMax 90." << std::endl ;
0264   std::cout <<  std::endl ;
0265   std::cout << "# use pseudo rapidity rather than polar angle - specify minimum and maximum eta value" << std::endl ;
0266   std::cout << "# etaMin -3." << std::endl ;
0267   std::cout << "# etaMax 3." << std::endl ;
0268   std::cout <<  std::endl ;
0269   std::cout << "# phi direction in deg (default: 90./y-axis)" << std::endl ;
0270   std::cout << "phi 90." << std::endl ;
0271   std::cout <<  std::endl ;
0272   std::cout << "# names and subdetector ranges given in [rmin,zmin,rmax,zmax] - e.g. for ILD_l5_vo2  (run dumpdetector -d to get numbers... ) " << std::endl ;
0273   std::cout <<  std::endl ;
0274   std::cout << "subdet vxd    0. 0. 6.549392e+00 1.450000e+01" << std::endl ;
0275   std::cout << "subdet sitftd 0. 0. 3.024000e+01 2.211800e+02" << std::endl ;
0276   std::cout << "subdet tpc    0. 0. 1.692100e+02 2.225000e+02" << std::endl ;
0277   std::cout << "subdet outtpc 0. 0. 1.769800e+02 2.350000e+02" << std::endl ;
0278   std::cout << "subdet set    0. 0. 1.775200e+02 2.300000e+02" << std::endl ;
0279   ::exit(0);
0280 }