GEAR  1.9.0
 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_terminate( gear_unexpected ) ;
55 
56  if( argc < 2 ) {
57  std::cout << " testgear: Testprogram for gear classes. " << std::endl
58  << " usage: testgear input.xml " << std::endl ;
59  exit(1) ;
60  }
61 
62  std::string fileName( argv[1] ) ;
63 
64  GearXML gearXML( fileName ) ;
65 
66  GearMgr* gearMgr = gearXML.createGearMgr() ;
67 
68  //convert all lengths to cm for TGeo
69 
70  const TPCParameters &tpcparams= gearMgr->getTPCParameters() ;
71  double r_in = tpcparams.getDoubleVal("tpcInnerRadius")/10.;
72  double r_o = tpcparams.getDoubleVal("tpcOuterRadius")/10.;
73  double wall_o = tpcparams.getDoubleVal("tpcOuterWallThickness")/10.;
74 
75  double Rtpc_i = r_in ; // -wall_in;
76  double Rtpc_o = r_o ; // + wall_o ;
77  double Rtpc_o_inner = r_o - wall_o ;
78  double Ztpc_inner = tpcparams.getMaxDriftLength()/10. ;
79 
80  const CalorimeterParameters &ecalparamsB= gearMgr->getEcalBarrelParameters() ;
81  double Recal_i = ecalparamsB.getExtent()[0]/10. - 3.0 ; // workaround for Ecal Si layer
82  // double Recal_i = ecalparamsB.getExtent()[0]/10. ; // workaround for Ecal Si layer
83 
84  const CalorimeterParameters &ecalparamsEC= gearMgr->getEcalEndcapParameters() ;
85  double Zecal=ecalparamsEC.getExtent()[2]/10;
86 
87  // Zecal = 235.1 ; // FIXME: exclude colling pipes for test !!!!!!
88  // Zecal = 246.0 ; // FIXME: include first ecal layers
89 
90  const ZPlanarParameters& vxdp = gearMgr->getVXDParameters() ;
91  double ro_VXD = vxdp.getShellOuterRadius() / 10. ;
92  double zo_VXD = vxdp.getShellHalfLength()/10. + vxdp.getShellGap()/10. ;
93 
95 
96  //===========================================================
97 
98  // for(int iTheta=0; iTheta<90; iTheta++ ) {
99  // double theta = iTheta * M_PI/180. ;
100  // double eta = - log( tan( theta / 2. ) ) ;
101 
102  // int N = 1000 ;
103  // double etaMin = 0. ; // 90 deg
104  // double etaMax = 3. ; // ~6 deg.
105  // double dEta = ( etaMax - etaMin ) / N ;
106 
107  // for(int i= 0 ; i < N ; ++i ) {
108  // double eta = etaMin + i * dEta ;
109  // double theta = 2. * atan( exp( -eta ) ) ;
110  // // eta = - log( tan( theta / 2. ) ) ;
111  // double iTheta = theta / M_PI * 180 ;
112 
113  int N = 10000 ;
114 
115  double thetaMin = 0. ;
116  double thetaMax = M_PI/2. ;
117 
118  // double thetaMin = 0.54 ;
119  // double thetaMax = 0.56 ;
120 
121  double dTheta = ( thetaMax - thetaMin ) / N ;
122 
123 
124  //std::cout << " iTheta, theta, eta, nrVXD, nrSIT, nrTPC, nrSET " << std::endl ;
125 
126  std::cerr << " ************** VXD cylinder r = " << ro_VXD << ", z = " << zo_VXD << std::endl ;
127  std::cerr << " ************** SIT cylinder r = " << Rtpc_i << ", z = " << Zecal << std::endl ;
128  std::cerr << " ************** TPC inner cylinder r = " << Rtpc_o_inner << ", z = " << Ztpc_inner << std::endl ;
129  std::cerr << " ************** TPC cylinder r = " << Rtpc_o << ", z = " << Zecal << std::endl ;
130  std::cerr << " ************** SET cylinder r = " << Recal_i << ", z = " << Zecal << std::endl ;
131 
132 
133  // write a String, so that you can directly read with root:
134  // TTree t ; t.ReadFile("intLen_ILD_o2_v05.txt")
135  std::cout << "itheta/d:theta/d:eta/d:nrvxd:nrsit/d:nrtpc_i/d:nrtpc/d:nrset/d:nrecal/d" << std::endl ;
136 
137  for(int i= 0 ; i < N ; ++i ) {
138  double theta = thetaMin + i * dTheta ;
139 
140  // double dTheta = 0.001 ; // mrad
141  // for(int theta= 0. ; theta < thetaMax ; theta += dTheta ) {
142 
143  double eta = ( theta == 0.0 ? 0. : - log( tan( theta / 2. ) ) );
144  double iTheta = theta / M_PI * 180 ;
145 
146 
147  //----------
148  Vector3D ip( 0.,0.,0. ) ;
149 
150  double nrVXD = distProp.getNRadlen( ip , pointOnCylinder( theta , ro_VXD , zo_VXD ) ) ;
151 
152  double nrSIT = distProp.getNRadlen( ip , pointOnCylinder( theta , Rtpc_i , Zecal ) ) ;
153 
154  double nrTPC_i = distProp.getNRadlen( ip , pointOnCylinder( theta , Rtpc_o_inner , Ztpc_inner ) ) ;
155 
156  double nrTPC = distProp.getNRadlen( ip , pointOnCylinder( theta , Rtpc_o , Zecal ) ) ;
157 
158  double nrSET = distProp.getNRadlen( ip , pointOnCylinder( theta , Recal_i , Zecal ) ) ;
159 
160 
161  // special treatment of ecal: we want the point after the firat absorber layer - regardless of whether this
162  // is in the barrel or endcap, i.e. we can't use a simple cylinder....
163  // Note: the values are taken from printMaterials ILD_o1_v05
164  double ecalL1_r = ( 1.85020000e+03 ) / 10. ;
165  double ecalL1_z = ( 2.45717250e+03 ) / 10. ;
166 
167  gear::Vector3D ecalL1Point = pointOnCylinder( theta , ecalL1_r , ecalL1_z ) ;
168 
169  if( ecalL1Point.z() > 2.350000000e+03/10. ){
170 
171  ecalL1Point = Vector3D( ecalL1_z * tan( theta ) , __PHI0 , ecalL1_z , Vector3D::cylindrical ) ;
172 
173  }
174 
175  double nrECal = distProp.getNRadlen( ip , ecalL1Point ) ;
176 
177 
178 
179  if( i== N-10 ) {
180 
181 
182  gear::Vector3D dir = pointOnCylinder( theta , Recal_i , Zecal ) - ip ;
183  dir = ( 1. / dir.r() ) * dir ;
184 
185 
186  std::vector<std::string> names =
187  distProp.getMaterialNames( ip , pointOnCylinder( theta , Recal_i , Zecal ) ) ;
188 
189  std::vector<double> radl ;
190 
191  std::vector<double> thicks =
192  distProp.getMaterialThicknesses( ip , pointOnCylinder( theta , Recal_i , Zecal ) ) ;
193 
194  double lambda = 0. ;
195 
196  std::cerr << " ################################# materials between the two points : \n "
197  << ip << "\n" << pointOnCylinder( theta , Recal_i , Zecal )
198  << " theta " << theta * 180. /3.141592 << "\n"
199  << " dir vector " << dir << "\n"
200  << " ############################################################# "
201  << std::endl ;
202 
203 
204 
205  for(unsigned j=0 ; j < names.size() ; ++j){
206 
207  lambda += thicks[j] ;
208 
209  gear::Vector3D local = ip + lambda * dir ;
210  double rL = distProp.getNRadlen( ip , local ) ;
211 
212  std::cerr << " ------ " << names[j] << " - " << thicks[j] << " cm - "
213  << rL
214  << std::endl ;
215 
216  }
217 
218 
219  }
220 
221  std::cout << iTheta << " "
222  << theta << " "
223  << eta << " "
224  << nrVXD << " "
225  << nrSIT << " "
226  << nrTPC_i << " "
227  << nrTPC << " "
228  << nrSET << " "
229  << nrECal
230  << std::endl;
231 
232  }
233 
234 
235 }
236 
237 
238 Vector3D pointOnCylinder( double theta, double r, double z, double phi){
239 
240  double theta0 = atan2( r, z ) ;
241 
242  // return ( theta > theta0 ?
243  // Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
244  // Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
245 
246  Vector3D v = ( theta > theta0 ?
247  Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
248  Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
249 
250 
251  // std::cout << " ==== pointOnCylinder( " << theta << ", " << r << ", " << z << ") " << " tan(theta) : " << tan( theta ) << std::endl
252  // << " " << v << std::endl ;
253 
254  return v ;
255 }
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.