Newer
Older
//==========================================================================
// 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
//
//==========================================================================
#include "TError.h"
// Framework include files
#include "DD4hep/Detector.h"
#include "DD4hep/DetType.h"
#include "DD4hep/Printout.h"
#include "DDRec/MaterialManager.h"
// #include "TGeoVolume.h"
// #include "TGeoManager.h"
// #include "TGeoNode.h"
#include "TFile.h"
#include "TH1F.h"
#include <fstream>
#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 };
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
};
/// comput a point on a cylinder (r,z) for a given direction (theta,phi) from the IP
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. ;
std::string outFileName("material_budget.root") ;
std::vector<SDetHelper> subdets ;
std::ifstream infile(steeringFile );
std::string 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" ){
iss >> det.name >> det.r0 >> det.z0 >> det.r1 >> det.z1 ;
subdets.emplace_back( det );
}
if ( !iss.eof() || iss.fail() ){
std::cout << " ERROR parsing line : " << line << std::endl ;
exit(1) ;
}
}
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 ;
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::cout << 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.),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 ) ;
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
std::cout << std::scientific << sum_x0 << " " << sum_lambda << " " ; // << path_length ;
}
std::cout << 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);
}