materialBudget.cpp

Go to the documentation of this file.

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
//==========================================================================
//
//  Compute the  material budget in terms of integrated radiation and ineraction
//  lengths as seen from the IP for various subdetector regions given in
//  a simple steering file. Creates a root file w/ histograms and dumps numbers
//  the screen.
// 
//  Author     : F.Gaede, DESY
//
//==========================================================================
 
// Framework include files
#include <DD4hep/DetType.h>
#include <DD4hep/Printout.h>
#include <DD4hep/Detector.h>
#include <DDRec/MaterialManager.h>
 
#include <TFile.h>
#include <TH1F.h>
#include <TError.h>
 
#include <cerrno>
#include <fstream>
#include <climits>
 
#include "main.h"
 
using namespace dd4hep;
using namespace dd4hep::rec;
 
void dumpExampleSteering() ;
 
namespace {
 
  struct SDetHelper{
    std::string name { };
    double r0 { 0.0 };
    double r1 { 0.0 };
    double z0 { 0.0 };
    double z1 { 0.0 };
    TH1F*  hx { nullptr };
    TH1F*  hl { nullptr };
  };
 
  Vector3D pointOnCylinder( double theta, double r, double z, double phi){
    
    double  theta0 = atan2( r, z ) ;
 
    Vector3D v = (   theta > theta0      ?
                     Vector3D(       r           , phi , r / tan( theta ) , Vector3D::cylindrical )  :
                     Vector3D( z * tan( theta )  , phi ,      z           , Vector3D::cylindrical )  ) ;
    return v ;
  }
 
}
 
int main_wrapper(int argc, char** argv)   {
  
  struct Handler  {
    Handler() { SetErrorHandler(Handler::print); }
    
    static void print(int level, Bool_t abort, const char *location, const char *msg)  {
      if ( level > kInfo || abort ) ::printf("%s: %s\n", location, msg);
    }
    
    static void usage()  {
      std::cout << " usage: materialBudget compact.xml steering.txt" << std::endl 
                << "     -> create histograms with the material budget as seen from the IP within fixed ranges of (rmin, rmax, zmin. zmax)" << std::endl 
                << "        see example steering file for details ..." << std::endl 
                << "  -x   : dump example steering file " << std::endl
                << std::endl;
      exit(1);
    }
 
  } _handler;
 
 
  if( argc == 2 ){
    std::string dummy = argv[1] ;
    
    if( dummy == "-x" )
      dumpExampleSteering() ;
    else
      Handler::usage();
  }
  
  if( argc != 3 )
    Handler::usage();
  
 
  // ================== parse steering file ========================================================
  std::string compactFile  =  argv[1];
  std::string steeringFile =  argv[2];
  int nbins = 90 ;
  double phi0 = M_PI / 2. ;
  double thetaMin = 0. ;
  double thetaMax = 90. ;
  double etaMin = 0. ;
  double etaMax = -1. ;
  std::string outFileName("material_budget.root") ;
  std::vector<SDetHelper> subdets ;
    
 
  std::ifstream infile(steeringFile );
  std::string line;
 
  while (std::getline(infile, line))  {
    // ignore empty lines and comments
    if( line.empty() || line.find("#") == 0 )
      continue ;
 
    std::istringstream iss(line);
    std::string token ;
    iss >> token ;
    
    if( token == "nbins" ){
      iss >> nbins ;
    }
    else if( token == "thetaMin" ){
      iss >> thetaMin ;
    }
    else if( token == "thetaMax" ){
      iss >> thetaMax ;
    }
    else if( token == "etaMin" ){
      iss >> etaMin ;
    }
    else if( token == "etaMax" ){
      iss >> etaMax ;
    }
    else if( token == "rootfile" ){
      iss >> outFileName ;
    }
    else if( token == "phi" ){
      iss >> phi0 ;
      phi0 = phi0 / 180. * M_PI ;
    }
    else if( token == "subdet" ){
      SDetHelper det;
      iss >>  det.name >> det.r0 >> det.z0 >> det.r1 >> det.z1 ;
      subdets.emplace_back( det );
    }
 
    if ( iss.fail() ){
      std::cout << " ERROR parsing line : " << line << std::endl ;
      ::exit(EINVAL);
    }    
  }
 
  if ( nbins <= 0 )    {
    std::cout << "Invalid value for binning (nbins) " << nbins << std::endl;
    return EINVAL;
  }
 
  // =================================================================================================
  setPrintLevel(WARNING);
 
  Detector& description = Detector::getInstance();
  description.fromXML(compactFile ) ;
 
  //----- open root file and book histograms
  TFile* rootFile = new TFile(outFileName.c_str(),"recreate");
  for( auto& det : subdets )  {
    std::string hxn(det.name), hxnn(det.name) ;
    std::string hln(det.name), hlnn(det.name) ;
 
    if( etaMax > 0. )  {  // use eta
      hxn += "x0" ;
      hxnn += " integrated X0 vs eta" ;
      det.hx = new TH1F( hxn.c_str(), hxnn.c_str(), nbins, etaMin , etaMax ) ;
 
      hln += "lambda" ;
      hlnn += " integrated int. lengths vs eta" ;
      det.hl = new TH1F( hln.c_str(), hlnn.c_str(), nbins, etaMin , etaMax ) ;
    }
    else {   // use polar angle
      hxn += "x0" ;
      hxnn += " integrated X0 vs -theta" ;
      det.hx = new TH1F( hxn.c_str(), hxnn.c_str(), nbins, -thetaMax , -thetaMin ) ;
 
      hln += "lambda" ;
      hlnn += " integrated int. lengths vs -theta" ;
      det.hl = new TH1F( hln.c_str(), hlnn.c_str(), nbins, -thetaMax , -thetaMin ) ;
    }
  }
  //-------------------------
      
 
  Volume world  = description.world().volume() ;
  
  MaterialManager matMgr( world ) ;
 
  thetaMin = thetaMin / 180. * M_PI ;
  thetaMax = thetaMax / 180. * M_PI ;
  double dTheta = (thetaMax-thetaMin)/nbins; // bin size
  double dEta  = (etaMax-etaMin)/nbins ;
 
  std::cout  << "====================================================================================================" << std::endl ;
  std::cout  << "theta:f/" ;
  for(auto& det : subdets){ std::cout  << det.name << "_x0:f/" << det.name << "_lam:f/" ; }
  std::cout  << std::endl ;
 
  if ( nbins <= 0 || nbins > USHRT_MAX )  {
    std::cout << "Unreasonable number of bins: " << nbins << std::endl;
    ::exit(EINVAL);
  }
 
  for(int i=0 ; i< nbins ;++i){
    double theta = ( etaMax > 0. ?  2. * atan ( exp ( - ( etaMin + (0.5+i)*dEta) ) ) : ( thetaMin + (0.5+i)*dTheta ) ) ;
    std::stringstream paramLine;
 
    paramLine << std::scientific << theta << " " ;
    for( auto& det : subdets )  {
      Vector3D p0 = pointOnCylinder( theta, det.r0 , det.z0 , phi0  ) ;// double theta, double r, double z, double phi)
      Vector3D p1 = pointOnCylinder( theta, det.r1 , det.z1 , phi0  ) ;// double theta, double r, double z, double phi)
      const MaterialVec& materials = matMgr.materialsBetween(p0, p1);
      double sum_x0(0.), sum_lambda(0.);
      // double path_length(0.);
 
      for( auto amat : materials )  {
        TGeoMaterial* mat =  amat.first->GetMaterial();
        double length = amat.second;
        double nx0 = length / mat->GetRadLen();
        sum_x0 += nx0;
        double nLambda = length / mat->GetIntLen();
        sum_lambda += nLambda;
        // path_length += length;
      }
 
      double binX = ( etaMax > 0. ? (etaMin + (0.5+i)*dEta) : -theta/M_PI*180. ) ;
      det.hx->Fill( binX , sum_x0 ) ;
      det.hl->Fill( binX , sum_lambda ) ;
      paramLine  << std::scientific  << sum_x0 << "  " << sum_lambda << "  " ; // << path_length ;
    }
    std::cout << paramLine.str() << std::endl;
  }  
  std::cout  << "====================================================================================================" << std::endl ;
  rootFile->Write();
  rootFile->Close();
  return 0;
}
 
void dumpExampleSteering(){
  std::cout << "# Example steering file for materialBudget  (taken from ILD_l5_v02)" << std::endl ;
  std::cout <<  std::endl ;
  std::cout << "# root output file"  << std::endl ;
  std::cout << "rootfile material_budget.root" << std::endl ;
  std::cout << "# number of bins for polar angle (default 90)" << std::endl ;
  std::cout << "nbins 90" << std::endl ;
  std::cout <<  std::endl ;
  std::cout << "# use polar angle - specify minimum and maximum theta value" << std::endl ;
  std::cout << "# thetaMin 0." << std::endl ;
  std::cout << "# thetaMax 90." << std::endl ;
  std::cout <<  std::endl ;
  std::cout << "# use pseudo rapidity rather than polar angle - specify minimum and maximum eta value" << std::endl ;
  std::cout << "# etaMin -3." << std::endl ;
  std::cout << "# etaMax 3." << std::endl ;
  std::cout <<  std::endl ;
  std::cout << "# phi direction in deg (default: 90./y-axis)" << std::endl ;
  std::cout << "phi 90." << std::endl ;
  std::cout <<  std::endl ;
  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 ;
  std::cout <<  std::endl ;
  std::cout << "subdet vxd    0. 0. 6.549392e+00 1.450000e+01" << std::endl ;
  std::cout << "subdet sitftd 0. 0. 3.024000e+01 2.211800e+02" << std::endl ;
  std::cout << "subdet tpc    0. 0. 1.692100e+02 2.225000e+02" << std::endl ;
  std::cout << "subdet outtpc 0. 0. 1.769800e+02 2.350000e+02" << std::endl ;
  std::cout << "subdet set    0. 0. 1.775200e+02 2.300000e+02" << std::endl ;
  ::exit(0);
}