GEAR  1.9.0
 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_terminate( gear_unexpected ) ;
104 
105  if( argc < 2 ) {
106  std::cout << " testgear: Testprogram for gear classes. " << std::endl
107  << " usage: testgear input.xml " << std::endl ;
108  exit(1) ;
109  }
110 
111  std::string fileName( argv[1] ) ;
112 
113  GearXML gearXML( fileName ) ;
114 
115  GearMgr* gearMgr = gearXML.createGearMgr() ;
116 
117  //convert all lengths to cm for TGeo
118 
119  // // const TPCParameters &tpcparams= gearMgr->getTPCParameters() ;
120  // // double r_in = tpcparams.getDoubleVal("tpcInnerRadius")/10.;
121  // // double r_o = tpcparams.getDoubleVal("tpcOuterRadius")/10.;
122  // // // double wall_o = tpcparams.getDoubleVal("tpcOuterWallThickness")/10.;
123 
124  // // double Rtpc_i = r_in ; // -wall_in;
125  // // double Rtpc_o = r_o ; // + wall_o ;
126 
127 
128  // float Cos8 = cos(M_PI/8.0);
129  float Cos12 = cos(M_PI/12.0);
130  float Cos16 = cos(M_PI/16.);
131 
132  const CalorimeterParameters &ecalparamsB= gearMgr->getEcalBarrelParameters() ;
133  double Recal_i = ecalparamsB.getExtent()[0]/10. ;
134  double Recal_o = ecalparamsB.getExtent()[1]/10. ;
135  //Recal_o /= Cos8 ; // we are at very small angles
136 
137  const CalorimeterParameters &ecalparamsEC= gearMgr->getEcalEndcapParameters() ;
138  double Zecal_i=ecalparamsEC.getExtent()[2]/10;
139  double Zecal_o=ecalparamsEC.getExtent()[3]/10;
140 
141  const CalorimeterParameters &hcalparamsB= gearMgr->getHcalBarrelParameters() ;
142  double Rhcal_i = hcalparamsB.getExtent()[0]/10. ;
143  double Rhcal_o = hcalparamsB.getExtent()[1]/10. ;
144  Rhcal_o /= Cos16 ;
145  const CalorimeterParameters &hcalparamsEC= gearMgr->getHcalEndcapParameters() ;
146  double Zhcal_i=hcalparamsEC.getExtent()[2]/10;
147  double Zhcal_o=hcalparamsEC.getExtent()[3]/10;
148 
149  const gear::GearParameters& parCoil = gearMgr->getGearParameters("CoilParameters");
150  double Zcoil_o = parCoil.getDoubleVal("Coil_cryostat_inner_cyl_half_z" )/10 ;
151  double Rcoil_i = parCoil.getDoubleVal("Coil_cryostat_inner_cyl_inner_radius" )/10 ;
152  double Rcoil_o = parCoil.getDoubleVal("Coil_cryostat_inner_cyl_outer_radius" )/10 ;
153 
154  const CalorimeterParameters &yokeparamsB= gearMgr->getYokeBarrelParameters() ;
155  double Ryoke_i = yokeparamsB.getExtent()[0]/10. ;
156  double Ryoke_o = yokeparamsB.getExtent()[1]/10. ;
157  Ryoke_o /= Cos12 ;
158  const CalorimeterParameters &yokeparamsEC= gearMgr->getYokeEndcapParameters() ;
159  double Zyoke_o=yokeparamsEC.getExtent()[3]/10;
160  const gear::CalorimeterParameters& yokeplug = gearMgr->getYokePlugParameters();
161  double Zyoke_i=yokeplug.getExtent()[2]/10;
162 
163 
164  // average between consecutive detectors (this should be right in the gaps)
165  Recal_i -= 20. ;
166  Zecal_i -= 20. ;
167 
168  Recal_o = ( Recal_o + Rhcal_i ) / 2. ;
169  Zecal_o = ( Zecal_o + Zhcal_i ) / 2. ;
170 
171  Rhcal_o = ( Rhcal_o + Rcoil_i ) / 2. ;
172  Zhcal_o = ( Zhcal_o + Zyoke_i ) / 2. ;
173 
174  Rcoil_o = ( Rcoil_o + Ryoke_i ) / 2. ;
175  Zcoil_o = Zhcal_o ;
176 
177  Ryoke_o += 50. ;
178  Zyoke_o += 50. ;
179 
180 
181 
183 
184  //===========================================================
185 
186  // for(int iTheta=0; iTheta<90; iTheta++ ) {
187  // double theta = iTheta * M_PI/180. ;
188  // double eta = - log( tan( theta / 2. ) ) ;
189 
190  // int N = 1000 ;
191  // double etaMin = 0. ; // 90 deg
192  // double etaMax = 3. ; // ~6 deg.
193  // double dEta = ( etaMax - etaMin ) / N ;
194 
195  // for(int i= 0 ; i < N ; ++i ) {
196  // double eta = etaMin + i * dEta ;
197  // double theta = 2. * atan( exp( -eta ) ) ;
198  // // eta = - log( tan( theta / 2. ) ) ;
199  // double iTheta = theta / M_PI * 180 ;
200 
201  int N = 1000 ;
202 
203  double thetaMin = 0. ;
204  double thetaMax = M_PI/2. ;
205 
206  // double thetaMin = 0.54 ;
207  // double thetaMax = 0.56 ;
208 
209  double dTheta = ( thetaMax - thetaMin ) / N ;
210 
211 
212  //std::cout << " iTheta, theta, eta, nrVXD, nrSIT, nrTPC, nrSET " << std::endl ;
213 
214 
215  std::cerr << " ************** SET cylinder r = " << Recal_i << ", z = " << Zecal_i << std::endl ;
216  std::cerr << " ************** Ecal cylinder r = " << Recal_o << ", z = " << Zecal_o << std::endl ;
217  std::cerr << " ************** Hcal cylinder r = " << Rhcal_o << ", z = " << Zhcal_o << std::endl ;
218  std::cerr << " ************** Coil cylinder r = " << Rcoil_o << ", z = " << Zcoil_o << std::endl ;
219  std::cerr << " ************** Yoke cylinder r = " << Ryoke_o << ", z = " << Zyoke_o << std::endl ;
220 
221 
222  // write a String, so that you can directly read with root:
223  // TTree t ; t.ReadFile("intLen_ILD_o2_v05.txt")
224 
225  std::cout << "iTheta/d:theta/d:eta/d:nilSET/d:nilEcal/d:nilHcal/d:nilCoil/d:nilYoke/d" << std::endl ;
226 
227  for(int i= 0 ; i < N ; ++i ) {
228  double theta = thetaMin + i * dTheta ;
229 
230  // double dTheta = 0.001 ; // mrad
231  // for(int theta= 0. ; theta < thetaMax ; theta += dTheta ) {
232 
233  double eta = ( theta == 0.0 ? 0. : - log( tan( theta / 2. ) ) );
234  double iTheta = theta / M_PI * 180 ;
235 
236  //----------
237  Vector3D ip( 0.,0.,0. ) ;
238 
239  Vector3D pSET = pointOnCylinder( theta , Recal_i , Zecal_i ) ;
240  Vector3D pEcal = pointOnCylinder( theta , Recal_o , Zecal_o ) ;
241  Vector3D pHcal = pointOnCylinder( theta , Rhcal_o , Zhcal_o ) ;
242  Vector3D pCoil = pointOnCylinder( theta , Rcoil_o, Zcoil_o ) ;
243  Vector3D pYoke = pointOnCylinder( theta , Ryoke_o , Zyoke_o ) ;
244 
245  double nilSET = distProp.getNIntlen( ip , pSET ) ;
246  double nilEcal = distProp.getNIntlen( ip , pEcal ) ;
247  double nilHcal = distProp.getNIntlen( ip , pHcal ) ;
248  double nilCoil = distProp.getNIntlen( ip , pCoil ) ;
249  double nilYoke = distProp.getNIntlen( ip , pYoke ) ;
250 
251  if( i == N-5 ) {
252 
253  printMaterial( distProp , pSET , pEcal ) ;
254  printMaterial( distProp , pEcal , pHcal ) ;
255  printMaterial( distProp , pHcal , pCoil ) ;
256  printMaterial( distProp , pCoil , pYoke ) ;
257 
258  // double R_o = Ryoke_i ;
259  // double Z_o = Zyoke_
260  // gear::Vector3D start = pointOnCylinder( theta , Recal_i , Zecal_i ) ;
261  // gear::Vector3D point = pointOnCylinder( theta , R_o , Z_o ) ;
262  // gear::Vector3D dir = point - start ;
263  // dir = ( 1. / dir.r() ) * dir ;
264  // std::vector<std::string> names =
265  // distProp.getMaterialNames( start , point ) ;
266  // std::vector<double> radl ;
267  // std::vector<double> thicks =
268  // distProp.getMaterialThicknesses( start , point ) ;
269  // double lambda = 0. ;
270  //
271  // std::cerr << " ################################# materials between the two points : \n "
272  // << start << "\n" << point
273  // << " theta " << theta * 180. /3.141592 << "\n"
274  // << " dir vector " << dir << "\n"
275  // << " ############################################################# "
276  // << std::endl ;
277  // for(unsigned j=0 ; j < names.size() ; ++j){
278  // lambda += thicks[j] ;
279  // gear::Vector3D local = start + lambda * dir ;
280  // double rL = distProp.getNRadlen( start, local ) ;
281  // double iL = distProp.getNIntlen( start, local ) ;
282  // std::cerr << " ------ " << names[j] << " r: " << lambda << " - " << thicks[j] << " cm - " << rL << " - iL: " << iL
283  // << std::endl ;
284  // }
285 
286  }
287 
288 
289  // ###################################################
290 
291  std::cout << iTheta << " "
292  << theta << " "
293  << eta << " "
294  << nilSET << " "
295  << nilEcal << " "
296  << nilHcal << " "
297  << nilCoil << " "
298  << nilYoke << " "
299  << std::endl;
300 
301 
302  // ###################################################
303 
304  }
305 
306 }
307 
308 
309 Vector3D pointOnCylinder( double theta, double r, double z, double phi){
310 
311  double theta0 = atan2( r, z ) ;
312 
313  // return ( theta > theta0 ?
314  // Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
315  // Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
316 
317  Vector3D v = ( theta > theta0 ?
318  Vector3D( r , phi , r / tan( theta ) , Vector3D::cylindrical ) :
319  Vector3D( z * tan( theta ) , phi , z , Vector3D::cylindrical ) ) ;
320 
321 
322  // std::cout << " ==== pointOnCylinder( " << theta << ", " << r << ", " << z << ") " << " tan(theta) : " << tan( theta ) << std::endl
323  // << " " << v << std::endl ;
324 
325  return v ;
326 }
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.