GEAR  1.6.1
 All Classes Namespaces Functions Variables Typedefs Enumerations Friends Pages
testMaterialBudgetNew.cc
1 
2 #include "gearimpl/Util.h"
3 #include "gearxml/GearXML.h"
4 #include "gear/GearMgr.h"
5 #include "gear/GEAR.h"
6 
7 #include "geartgeo/TGeoGearDistanceProperties.h"
8 #include "geartgeo/TGeoGearPointProperties.h"
9 
10 #include <iostream>
11 #include <assert.h>
12 
13 #include <exception>
14 #include <typeinfo>
15 #include <cstdlib>
16 
17 #include <sstream>
18 #include <fstream>
19 #include <math.h>
20 
21 using namespace gear ;
22 
23 #define __PHI0 0.0
24 
25 void gear_unexpected(){
26 
27  try {
28 
29  throw ;
30 
31  } catch( std::exception& e) {
32 
33  std::cout << " A runtime error has occured : "
34  << e.what()
35  << std::endl
36  << " the program will have to be terminated - sorry." << std::endl ;
37  exit(1) ;
38  }
39 }
40 
41 
42 
43 Vector3D pointOnCylinder( double theta, double r, double z, double phi= __PHI0 ) ;
44 
45 
51 int main(int argc, char**argv){
52 
53 
54  std::set_unexpected( gear_unexpected ) ;
55  std::set_terminate( gear_unexpected ) ;
56 
57  if( argc < 2 ) {
58  std::cout << " testgear: Testprogram for gear classes. " << std::endl
59  << " usage: testgear input.xml " << std::endl ;
60  exit(1) ;
61  }
62 
63  std::string fileName( argv[1] ) ;
64 
65  GearXML gearXML( fileName ) ;
66 
67  GearMgr* gearMgr = gearXML.createGearMgr() ;
68 
69  //convert all lengths to cm for TGeo
70 
71  const TPCParameters &tpcparams= gearMgr->getTPCParameters() ;
72  double r_in = tpcparams.getDoubleVal("tpcInnerRadius")/10.;
73  double r_o = tpcparams.getDoubleVal("tpcOuterRadius")/10.;
74  double wall_o = tpcparams.getDoubleVal("tpcOuterWallThickness")/10.;
75 
76  double Rtpc_i = r_in ; // -wall_in;
77  double Rtpc_o = r_o ; // + wall_o ;
78  double Rtpc_o_inner = r_o - wall_o ;
79  double Ztpc_inner = tpcparams.getMaxDriftLength()/10. ;
80 
81  const CalorimeterParameters &ecalparamsB= gearMgr->getEcalBarrelParameters() ;
82  double Recal_i = ecalparamsB.getExtent()[0]/10. - 3.0 ; // workaround for Ecal Si layer
83  // double Recal_i = ecalparamsB.getExtent()[0]/10. ; // workaround for Ecal Si layer
84 
85  const CalorimeterParameters &ecalparamsEC= gearMgr->getEcalEndcapParameters() ;
86  double Zecal=ecalparamsEC.getExtent()[2]/10;
87 
88  // Zecal = 235.1 ; // FIXME: exclude colling pipes for test !!!!!!
89  // Zecal = 246.0 ; // FIXME: include first ecal layers
90 
91  const ZPlanarParameters& vxdp = gearMgr->getVXDParameters() ;
92  double ro_VXD = vxdp.getShellOuterRadius() / 10. ;
93  double zo_VXD = vxdp.getShellHalfLength()/10. + vxdp.getShellGap()/10. ;
94 
96 
97  //===========================================================
98 
99  // for(int iTheta=0; iTheta<90; iTheta++ ) {
100  // double theta = iTheta * M_PI/180. ;
101  // double eta = - log( tan( theta / 2. ) ) ;
102 
103  // int N = 1000 ;
104  // double etaMin = 0. ; // 90 deg
105  // double etaMax = 3. ; // ~6 deg.
106  // double dEta = ( etaMax - etaMin ) / N ;
107 
108  // for(int i= 0 ; i < N ; ++i ) {
109  // double eta = etaMin + i * dEta ;
110  // double theta = 2. * atan( exp( -eta ) ) ;
111  // // eta = - log( tan( theta / 2. ) ) ;
112  // double iTheta = theta / M_PI * 180 ;
113 
114  int N = 10000 ;
115 
116  double thetaMin = 0. ;
117  double thetaMax = M_PI/2. ;
118 
119  // double thetaMin = 0.54 ;
120  // double thetaMax = 0.56 ;
121 
122  double dTheta = ( thetaMax - thetaMin ) / N ;
123 
124 
125  //std::cout << " iTheta, theta, eta, nrVXD, nrSIT, nrTPC, nrSET " << std::endl ;
126 
127  std::cerr << " ************** VXD cylinder r = " << ro_VXD << ", z = " << zo_VXD << std::endl ;
128  std::cerr << " ************** SIT cylinder r = " << Rtpc_i << ", z = " << Zecal << std::endl ;
129  std::cerr << " ************** TPC inner cylinder r = " << Rtpc_o_inner << ", z = " << Ztpc_inner << std::endl ;
130  std::cerr << " ************** TPC cylinder r = " << Rtpc_o << ", z = " << Zecal << std::endl ;
131  std::cerr << " ************** SET cylinder r = " << Recal_i << ", z = " << Zecal << std::endl ;
132 
133 
134  // write a String, so that you can directly read with root:
135  // TTree t ; t.ReadFile("intLen_ILD_o2_v05.txt")
136  std::cout << "itheta/d:theta/d:eta/d:nrvxd:nrsit/d:nrtpc_i/d:nrtpc/d:nrset/d:nrecal/d" << std::endl ;
137 
138  for(int i= 0 ; i < N ; ++i ) {
139  double theta = thetaMin + i * dTheta ;
140 
141  // double dTheta = 0.001 ; // mrad
142  // for(int theta= 0. ; theta < thetaMax ; theta += dTheta ) {
143 
144  double eta = ( theta == 0.0 ? 0. : - log( tan( theta / 2. ) ) );
145  double iTheta = theta / M_PI * 180 ;
146 
147 
148  //----------
149  Vector3D ip( 0.,0.,0. ) ;
150 
151  double nrVXD = distProp.getNRadlen( ip , pointOnCylinder( theta , ro_VXD , zo_VXD ) ) ;
152 
153  double nrSIT = distProp.getNRadlen( ip , pointOnCylinder( theta , Rtpc_i , Zecal ) ) ;
154 
155  double nrTPC_i = distProp.getNRadlen( ip , pointOnCylinder( theta , Rtpc_o_inner , Ztpc_inner ) ) ;
156 
157  double nrTPC = distProp.getNRadlen( ip , pointOnCylinder( theta , Rtpc_o , Zecal ) ) ;
158 
159  double nrSET = distProp.getNRadlen( ip , pointOnCylinder( theta , Recal_i , Zecal ) ) ;
160 
161 
162  // special treatment of ecal: we want the point after the firat absorber layer - regardless of whether this
163  // is in the barrel or endcap, i.e. we can't use a simple cylinder....
164  // Note: the values are taken from printMaterials ILD_o1_v05
165  double ecalL1_r = ( 1.85020000e+03 ) / 10. ;
166  double ecalL1_z = ( 2.45717250e+03 ) / 10. ;
167 
168  gear::Vector3D ecalL1Point = pointOnCylinder( theta , ecalL1_r , ecalL1_z ) ;
169 
170  if( ecalL1Point.z() > 2.350000000e+03/10. ){
171 
172  ecalL1Point = Vector3D( ecalL1_z * tan( theta ) , __PHI0 , ecalL1_z , Vector3D::cylindrical ) ;
173 
174  }
175 
176  double nrECal = distProp.getNRadlen( ip , ecalL1Point ) ;
177 
178 
179 
180  if( i== N-10 ) {
181 
182 
183  gear::Vector3D dir = pointOnCylinder( theta , Recal_i , Zecal ) - ip ;
184  dir = ( 1. / dir.r() ) * dir ;
185 
186 
187  std::vector<std::string> names =
188  distProp.getMaterialNames( ip , pointOnCylinder( theta , Recal_i , Zecal ) ) ;
189 
190  std::vector<double> radl ;
191 
192  std::vector<double> thicks =
193  distProp.getMaterialThicknesses( ip , pointOnCylinder( theta , Recal_i , Zecal ) ) ;
194 
195  double lambda = 0. ;
196 
197  std::cerr << " ################################# materials between the two points : \n "
198  << ip << "\n" << pointOnCylinder( theta , Recal_i , Zecal )
199  << " theta " << theta * 180. /3.141592 << "\n"
200  << " dir vector " << dir << "\n"
201  << " ############################################################# "
202  << std::endl ;
203 
204 
205 
206  for(unsigned j=0 ; j < names.size() ; ++j){
207 
208  lambda += thicks[j] ;
209 
210  gear::Vector3D local = ip + lambda * dir ;
211  double rL = distProp.getNRadlen( ip , local ) ;
212 
213  std::cerr << " ------ " << names[j] << " - " << thicks[j] << " cm - "
214  << rL
215  << std::endl ;
216 
217  }
218 
219 
220  }
221 
222  std::cout << iTheta << " "
223  << theta << " "
224  << eta << " "
225  << nrVXD << " "
226  << nrSIT << " "
227  << nrTPC_i << " "
228  << nrTPC << " "
229  << nrSET << " "
230  << nrECal
231  << std::endl;
232 
233  }
234 
235 
236 }
237 
238 
239 Vector3D pointOnCylinder( double theta, double r, double z, double phi){
240 
241  double theta0 = atan2( r, z ) ;
242 
243  // return ( theta > theta0 ?
244  // Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
245  // Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
246 
247  Vector3D v = ( theta > theta0 ?
248  Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
249  Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
250 
251 
252  // std::cout << " ==== pointOnCylinder( " << theta << ", " << r << ", " << z << ") " << " tan(theta) : " << tan( theta ) << std::endl
253  // << " " << v << std::endl ;
254 
255  return v ;
256 }
virtual double getShellOuterRadius() const =0
The outer radius of the support shell in mm.
Geometry properties of a vertex detector needed for reconstruction code.
virtual const TPCParameters & getTPCParameters() const =0
Get the TPCParameters.
Proposal for an abstract interface that defines the geometry properties of a TPC like detector needed...
Definition: TPCParameters.h:24
virtual const ZPlanarParameters & getVXDParameters() const =0
Get the VXD parameters.
virtual double getNRadlen(const Vector3D &p0, const Vector3D &p1) const
The number of radiation lengths along the distance between [p0,p1] .
virtual const std::vector< std::string > & getMaterialNames(const Vector3D &p0, const Vector3D &p1) const
List of matrial names along the distance between [p0,p1]- WARNING: this method returns a reference to...
virtual const GearDistanceProperties & getDistanceProperties() const =0
Get the distance properties object.
Proposal for an abstract interface that defines geometry properties of a typical sampling calorimeter...
TGeo Implementation of the abstract interface that returns the (material) properties along a given di...
virtual double getShellGap() const =0
The length of the gap in mm (gap position at z=0)
virtual const std::vector< double > & getExtent() const =0
Extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
Implementation of GEAR using XML.
Definition: GearXML.h:18
Simple three dimensional vector providing the components for cartesian, cylindrical and spherical coo...
Definition: Vector3D.h:18
virtual const std::vector< double > & getMaterialThicknesses(const Vector3D &p0, const Vector3D &p1) const
List of matrial thicknesses in mm along the distance between [p0,p1] - runs parallel to the array ret...
double z()
Cartesian cartesian z coordinate.
Definition: Vector3D.h:62
virtual const CalorimeterParameters & getEcalBarrelParameters() const =0
Get the Ecal barrel parameters.
virtual double getShellHalfLength() const =0
The half length (z) of the support shell in mm (w/o gap).
double r() const
Spherical r/magnitude.
Definition: Vector3D.h:119
Abstract interface for a manager class that returns the Gear classes for the relevant subdetectors...
Definition: GearMgr.h:36
virtual const CalorimeterParameters & getEcalEndcapParameters() const =0
Get the Ecal endcap parameters.
virtual double getMaxDriftLength() const =0
The maximum drift length in the TPC in mm.
virtual double getDoubleVal(const std::string &key) const =0
Double value for key.