GEAR  1.6.1
 All Classes Namespaces Functions Variables Typedefs Enumerations Friends Pages
testInteractionLengths.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 
24 void gear_unexpected(){
25 
26  try {
27 
28  throw ;
29 
30  } catch( std::exception& e) {
31 
32  std::cout << " A runtime error has occured : "
33  << e.what()
34  << std::endl
35  << " the program will have to be terminated - sorry." << std::endl ;
36  exit(1) ;
37  }
38 }
39 
40 
41 
42 Vector3D pointOnCylinder( double theta, double r, double z, double phi=0.02 ) ; //0.0) ;
43 
44 
45 void printMaterial( TGeoGearDistanceProperties& distProp, const gear::Vector3D& start , const gear::Vector3D& point ) {
46 
47 
48  gear::Vector3D dir = point - start ;
49  dir = ( 1. / dir.r() ) * dir ;
50 
51 
52  std::vector<std::string> names =
53  distProp.getMaterialNames( start , point ) ;
54 
55  std::vector<double> radl ;
56 
57  std::vector<double> thicks =
58  distProp.getMaterialThicknesses( start , point ) ;
59 
60  double lambda = 0. ;
61 
62  std::cerr << " ################################# materials between the two points : \n "
63  << start << "\n" << point
64  << " dir vector " << dir << "\n"
65  << " ############################################################# "
66  << std::endl ;
67 
68 
69 
70  for(unsigned j=0 ; j < names.size() ; ++j){
71 
72  lambda += thicks[j] ;
73 
74  gear::Vector3D local = start + lambda * dir ;
75  double rL = distProp.getNRadlen( start, local ) ;
76  double iL = distProp.getNIntlen( start, local ) ;
77 
78  std::cerr << " ------ " << names[j] << " r: " << lambda << " - " << thicks[j] << " cm - " << rL << " - iL: " << iL
79  << std::endl ;
80 
81  }
82 
83  std::cerr << "\n ================================================================================\n \n " ;
84 
85 }
86 
87 
88 
100 int main(int argc, char**argv){
101 
102 
103  std::set_unexpected( gear_unexpected ) ;
104  std::set_terminate( gear_unexpected ) ;
105 
106  if( argc < 2 ) {
107  std::cout << " testgear: Testprogram for gear classes. " << std::endl
108  << " usage: testgear input.xml " << std::endl ;
109  exit(1) ;
110  }
111 
112  std::string fileName( argv[1] ) ;
113 
114  GearXML gearXML( fileName ) ;
115 
116  GearMgr* gearMgr = gearXML.createGearMgr() ;
117 
118  //convert all lengths to cm for TGeo
119 
120  // // const TPCParameters &tpcparams= gearMgr->getTPCParameters() ;
121  // // double r_in = tpcparams.getDoubleVal("tpcInnerRadius")/10.;
122  // // double r_o = tpcparams.getDoubleVal("tpcOuterRadius")/10.;
123  // // // double wall_o = tpcparams.getDoubleVal("tpcOuterWallThickness")/10.;
124 
125  // // double Rtpc_i = r_in ; // -wall_in;
126  // // double Rtpc_o = r_o ; // + wall_o ;
127 
128 
129  float Cos8 = cos(M_PI/8.0);
130  float Cos12 = cos(M_PI/12.0);
131  float Cos16 = cos(M_PI/16.);
132 
133  const CalorimeterParameters &ecalparamsB= gearMgr->getEcalBarrelParameters() ;
134  double Recal_i = ecalparamsB.getExtent()[0]/10. ;
135  double Recal_o = ecalparamsB.getExtent()[1]/10. ;
136  //Recal_o /= Cos8 ; // we are at very small angles
137 
138  const CalorimeterParameters &ecalparamsEC= gearMgr->getEcalEndcapParameters() ;
139  double Zecal_i=ecalparamsEC.getExtent()[2]/10;
140  double Zecal_o=ecalparamsEC.getExtent()[3]/10;
141 
142  const CalorimeterParameters &hcalparamsB= gearMgr->getHcalBarrelParameters() ;
143  double Rhcal_i = hcalparamsB.getExtent()[0]/10. ;
144  double Rhcal_o = hcalparamsB.getExtent()[1]/10. ;
145  Rhcal_o /= Cos16 ;
146  const CalorimeterParameters &hcalparamsEC= gearMgr->getHcalEndcapParameters() ;
147  double Zhcal_i=hcalparamsEC.getExtent()[2]/10;
148  double Zhcal_o=hcalparamsEC.getExtent()[3]/10;
149 
150  const gear::GearParameters& pCoil = gearMgr->getGearParameters("CoilParameters");
151  double Zcoil_o = pCoil.getDoubleVal("Coil_cryostat_inner_cyl_half_z" )/10 ;
152  double Rcoil_i = pCoil.getDoubleVal("Coil_cryostat_inner_cyl_inner_radius" )/10 ;
153  double Rcoil_o = pCoil.getDoubleVal("Coil_cryostat_inner_cyl_outer_radius" )/10 ;
154 
155  const CalorimeterParameters &yokeparamsB= gearMgr->getYokeBarrelParameters() ;
156  double Ryoke_i = yokeparamsB.getExtent()[0]/10. ;
157  double Ryoke_o = yokeparamsB.getExtent()[1]/10. ;
158  Ryoke_o /= Cos12 ;
159  const CalorimeterParameters &yokeparamsEC= gearMgr->getYokeEndcapParameters() ;
160  double Zyoke_o=yokeparamsEC.getExtent()[3]/10;
161  const gear::CalorimeterParameters& yokeplug = gearMgr->getYokePlugParameters();
162  double Zyoke_i=yokeplug.getExtent()[2]/10;
163 
164 
165  // average between consecutive detectors (this should be right in the gaps)
166  Recal_i -= 20. ;
167  Zecal_i -= 20. ;
168 
169  Recal_o = ( Recal_o + Rhcal_i ) / 2. ;
170  Zecal_o = ( Zecal_o + Zhcal_i ) / 2. ;
171 
172  Rhcal_o = ( Rhcal_o + Rcoil_i ) / 2. ;
173  Zhcal_o = ( Zhcal_o + Zyoke_i ) / 2. ;
174 
175  Rcoil_o = ( Rcoil_o + Ryoke_i ) / 2. ;
176  Zcoil_o = Zhcal_o ;
177 
178  Ryoke_o += 50. ;
179  Zyoke_o += 50. ;
180 
181 
182 
184 
185  //===========================================================
186 
187  // for(int iTheta=0; iTheta<90; iTheta++ ) {
188  // double theta = iTheta * M_PI/180. ;
189  // double eta = - log( tan( theta / 2. ) ) ;
190 
191  // int N = 1000 ;
192  // double etaMin = 0. ; // 90 deg
193  // double etaMax = 3. ; // ~6 deg.
194  // double dEta = ( etaMax - etaMin ) / N ;
195 
196  // for(int i= 0 ; i < N ; ++i ) {
197  // double eta = etaMin + i * dEta ;
198  // double theta = 2. * atan( exp( -eta ) ) ;
199  // // eta = - log( tan( theta / 2. ) ) ;
200  // double iTheta = theta / M_PI * 180 ;
201 
202  int N = 1000 ;
203 
204  double thetaMin = 0. ;
205  double thetaMax = M_PI/2. ;
206 
207  // double thetaMin = 0.54 ;
208  // double thetaMax = 0.56 ;
209 
210  double dTheta = ( thetaMax - thetaMin ) / N ;
211 
212 
213  //std::cout << " iTheta, theta, eta, nrVXD, nrSIT, nrTPC, nrSET " << std::endl ;
214 
215 
216  std::cerr << " ************** SET cylinder r = " << Recal_i << ", z = " << Zecal_i << std::endl ;
217  std::cerr << " ************** Ecal cylinder r = " << Recal_o << ", z = " << Zecal_o << std::endl ;
218  std::cerr << " ************** Hcal cylinder r = " << Rhcal_o << ", z = " << Zhcal_o << std::endl ;
219  std::cerr << " ************** Coil cylinder r = " << Rcoil_o << ", z = " << Zcoil_o << std::endl ;
220  std::cerr << " ************** Yoke cylinder r = " << Ryoke_o << ", z = " << Zyoke_o << std::endl ;
221 
222 
223  // write a String, so that you can directly read with root:
224  // TTree t ; t.ReadFile("intLen_ILD_o2_v05.txt")
225 
226  std::cout << "iTheta/d:theta/d:eta/d:nilSET/d:nilEcal/d:nilHcal/d:nilCoil/d:nilYoke/d" << std::endl ;
227 
228  for(int i= 0 ; i < N ; ++i ) {
229  double theta = thetaMin + i * dTheta ;
230 
231  // double dTheta = 0.001 ; // mrad
232  // for(int theta= 0. ; theta < thetaMax ; theta += dTheta ) {
233 
234  double eta = ( theta == 0.0 ? 0. : - log( tan( theta / 2. ) ) );
235  double iTheta = theta / M_PI * 180 ;
236 
237  //----------
238  Vector3D ip( 0.,0.,0. ) ;
239 
240  Vector3D pSET = pointOnCylinder( theta , Recal_i , Zecal_i ) ;
241  Vector3D pEcal = pointOnCylinder( theta , Recal_o , Zecal_o ) ;
242  Vector3D pHcal = pointOnCylinder( theta , Rhcal_o , Zhcal_o ) ;
243  Vector3D pCoil = pointOnCylinder( theta , Rcoil_o, Zcoil_o ) ;
244  Vector3D pYoke = pointOnCylinder( theta , Ryoke_o , Zyoke_o ) ;
245 
246  double nilSET = distProp.getNIntlen( ip , pSET ) ;
247  double nilEcal = distProp.getNIntlen( ip , pEcal ) ;
248  double nilHcal = distProp.getNIntlen( ip , pHcal ) ;
249  double nilCoil = distProp.getNIntlen( ip , pCoil ) ;
250  double nilYoke = distProp.getNIntlen( ip , pYoke ) ;
251 
252  if( i == N-5 ) {
253 
254  printMaterial( distProp , pSET , pEcal ) ;
255  printMaterial( distProp , pEcal , pHcal ) ;
256  printMaterial( distProp , pHcal , pCoil ) ;
257  printMaterial( distProp , pCoil , pYoke ) ;
258 
259  // double R_o = Ryoke_i ;
260  // double Z_o = Zyoke_
261  // gear::Vector3D start = pointOnCylinder( theta , Recal_i , Zecal_i ) ;
262  // gear::Vector3D point = pointOnCylinder( theta , R_o , Z_o ) ;
263  // gear::Vector3D dir = point - start ;
264  // dir = ( 1. / dir.r() ) * dir ;
265  // std::vector<std::string> names =
266  // distProp.getMaterialNames( start , point ) ;
267  // std::vector<double> radl ;
268  // std::vector<double> thicks =
269  // distProp.getMaterialThicknesses( start , point ) ;
270  // double lambda = 0. ;
271  //
272  // std::cerr << " ################################# materials between the two points : \n "
273  // << start << "\n" << point
274  // << " theta " << theta * 180. /3.141592 << "\n"
275  // << " dir vector " << dir << "\n"
276  // << " ############################################################# "
277  // << std::endl ;
278  // for(unsigned j=0 ; j < names.size() ; ++j){
279  // lambda += thicks[j] ;
280  // gear::Vector3D local = start + lambda * dir ;
281  // double rL = distProp.getNRadlen( start, local ) ;
282  // double iL = distProp.getNIntlen( start, local ) ;
283  // std::cerr << " ------ " << names[j] << " r: " << lambda << " - " << thicks[j] << " cm - " << rL << " - iL: " << iL
284  // << std::endl ;
285  // }
286 
287  }
288 
289 
290  // ###################################################
291 
292  std::cout << iTheta << " "
293  << theta << " "
294  << eta << " "
295  << nilSET << " "
296  << nilEcal << " "
297  << nilHcal << " "
298  << nilCoil << " "
299  << nilYoke << " "
300  << std::endl;
301 
302 
303  // ###################################################
304 
305  }
306 
307 }
308 
309 
310 Vector3D pointOnCylinder( double theta, double r, double z, double phi){
311 
312  double theta0 = atan2( r, z ) ;
313 
314  // return ( theta > theta0 ?
315  // Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
316  // Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
317 
318  Vector3D v = ( theta > theta0 ?
319  Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
320  Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
321 
322 
323  // std::cout << " ==== pointOnCylinder( " << theta << ", " << r << ", " << z << ") " << " tan(theta) : " << tan( theta ) << std::endl
324  // << " " << v << std::endl ;
325 
326  return v ;
327 }
virtual const CalorimeterParameters & getHcalEndcapParameters() const =0
Get the Hcal endcap parameters.
virtual double getNRadlen(const Vector3D &p0, const Vector3D &p1) const
The number of radiation lengths along the distance between [p0,p1] .
virtual double getNIntlen(const Vector3D &p0, const Vector3D &p1) const
The number of interaction 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 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
virtual const CalorimeterParameters & getYokeBarrelParameters() const =0
Get the Yoke barrel parameters.
Abstract interface for a set of parameters that can be used to describe the geometrical properties of...
Simple three dimensional vector providing the components for cartesian, cylindrical and spherical coo...
Definition: Vector3D.h:18
virtual const GearParameters & getGearParameters(const std::string &key) const =0
Get named parameters for key.
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...
virtual const CalorimeterParameters & getEcalBarrelParameters() const =0
Get the Ecal barrel parameters.
virtual const CalorimeterParameters & getYokePlugParameters() const =0
Get the Yoke plug parameters.
virtual const CalorimeterParameters & getYokeEndcapParameters() const =0
Get the Yoke endcap parameters.
virtual const CalorimeterParameters & getHcalBarrelParameters() const =0
Get the Hcal barrel parameters.
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 getDoubleVal(const std::string &key) const =0
Double value for key.