GEAR  1.9.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Friends Pages
testVTXgear.cc
1 
2 #include "gearimpl/Util.h"
3 #include "gearxml/GearXML.h"
4 #include "gear/GearMgr.h"
5 #include "gear/GEAR.h"
6 #include "gear/VXDParameters.h"
7 
8 #ifdef CGA
9 #include "gearcga/CGAGearDistanceProperties.h"
10 #include "gearcga/CGAGearPointProperties.h"
11 #endif
12 
13 #include <iostream>
14 #include <assert.h>
15 
16 #include <exception>
17 #include <cstdlib>
18 
19 #include <sstream>
20 #include <fstream>
21 
22 #include <math.h>
23 
24 #ifdef GEAR_USE_AIDA
25 #include "AIDA/AIDA.h"
26 using namespace AIDA ;
27 #endif
28 
29 
30 
31 using namespace gear ;
32 
33 void testVXD( const VXDParameters& vxdParams ) ;
34 
35 void createCheckPlots( const VXDParameters& vxdParams ) ;
36 
37 void testVXDPoint( const Vector3D p , const VXDParameters& vxdParams ) ;
38 
39 void testVXDDist( const Vector3D p , const VXDParameters& vxdParams ) ;
40 
41 void testAllVXDPoint( const VXDParameters& vxdParams ) ;
42 
43 void testAllVXDDist( const VXDParameters& vxdParams ) ;
44 
45 void testVXDIntersection( const VXDParameters& vxdParams ) ;
46 
47 void gear_unexpected(){
48 
49  try {
50 
51  throw ;
52 
53  } catch( std::exception& e) {
54 
55  std::cout << " A runtime error has occured : "
56  << e.what()
57  << std::endl
58  << " the program will have to be terminated - sorry." << std::endl ;
59  exit(1) ;
60  }
61 }
62 
63 
64 #ifdef GEAR_USE_AIDA
65 // global variables for AIDA histograms and tuples
66 
67  IAnalysisFactory * myaida ;
68  ITreeFactory * mytreefactory ;
69  ITree * mytree ;
70  IHistogramFactory * myhistofactory ;
71  ITupleFactory * mytuplefactory ;
72 
73 #endif
74 
75 
76 
77 
78 
83 int main(int argc, char**argv){
84 
85 
86  std::set_terminate( gear_unexpected ) ;
87 
88  if( argc < 2 ) {
89  std::cout << " testgear: Testprogram for vertex detecto in gear. " << std::endl
90  << " usage: testgear input.xml " << std::endl ;
91  exit(1) ;
92  }
93 
94  std::string fileName( argv[1] ) ;
95 
96  GearXML gearXML( fileName ) ;
97 
98  GearMgr* gearMgr = gearXML.createGearMgr() ;
99 
100 
101  std::cout << " testVTXgear - instantiated GearMgr from file " << fileName
102  << std::endl ;
103 
104  const VXDParameters& vp =
105  dynamic_cast<const VXDParameters&>( gearMgr->getVXDParameters() ) ;
106 
107 
108  std::cout << " VXD parameters : " << std::endl
109  << vp
110  << std::endl ;
111 
112 
113 
114 #ifdef GEAR_USE_AIDA
115 
116  std::string storeName("vxdCheckPlot.root");
117 
118  myaida = AIDA_createAnalysisFactory() ;
119 
120  mytreefactory = myaida->createTreeFactory() ;
121 
122  mytree = mytreefactory->create(storeName,"root",false,true,"none") ;
123  // mytreefactory->create(storeName,storeType,readOnly,createNew,options)
124 
125  myhistofactory = myaida->createHistogramFactory( *mytree ) ;
126  mytuplefactory = myaida->createTupleFactory( *mytree ) ;
127 
128 // const VXDParameters& vp =
129 // dynamic_cast<const VXDParameters&>( gearMgr->getVXDParameters() ) ;
130 
131  //testVXDIntersection( vp ) ;
132 
133  createCheckPlots( vp ) ;
134 
135 #else
136 
137  std::cout << " Compile with GEAR_USE_AIDA==1 to get some histograms ! "
138  << std::endl;
139 
140 #endif
141 
142 
143  testVXD( vp ) ;
144 
145 }
146 
147 void testVXD( const VXDParameters& vxdParams ) {
148 
149  // Methods for special debugging
150 
151  Vector3D p0 ( 0, 0 , 0 ) ;
152  Vector3D p1 ( 12 , 2 , 0 ) ;
153 // Vector3D p1 ( -36.8 , 7 , 10 ) ;
154  Vector3D p2 ( -15 ,-8, .5 ) ;
155  Vector3D p3 ( -30 , -2 , .5 ) ;
156  Vector3D p4 ( -25 ,-25.0 , 0 ) ;
157 
158 // testVXDPoint( p0 , vxdParams ) ;
159 // testVXDPoint( p1 , vxdParams ) ;
160 // testVXDPoint( p2 , vxdParams ) ;
161 // testVXDPoint( p3 , vxdParams ) ;
162 // testVXDPoint( p4 , vxdParams ) ;
163 
164 
165  testVXDDist( p0, vxdParams ) ;
166  testVXDDist( p1, vxdParams ) ;
167  testVXDDist( p2, vxdParams ) ;
168  testVXDDist( p3, vxdParams ) ;
169  testVXDDist( p4, vxdParams ) ;
170 
171  // testAllVXDDist( vxdParams ) ;
172  testAllVXDPoint( vxdParams ) ;
173 }
174 
175 #ifdef GEAR_USE_AIDA
176 
177 void createCheckPlots( const VXDParameters& vxdParams ) {
178 
179  // creates sophisticated checkplots
180 
181  // values for LoRes Overview
182  double testStart = -65 ;
183  double testEnd = 65 ;
184  double testStep = .05 ;
185  double zStart = 10 ;
186  double zEnd = 11 ;
187  double zStep = 1 ;
188 
189  // Every _skipDist_ Step the distance test is performed
190  int skipDist = 1000 ;
191 
192  // values for HiRes Part
193  double testDetailStart = 2 ;
194  double testDetailEnd = 18 ;
195  double testDetailStep = 0.005 ;
196 
197  float myStatusP = 0, myStatusL = 0;
198  int nLadder = 0 , nSensitive = 0 ;
199 
200  // ladder and sensitive Map
201  std::string mapName = "vxdMap" ;
202  std::string mapTitle = "Map of the VTX" ;
203  std::vector<std::string> columnNames ;
204  std::vector<std::string> columnTypes ;
205 
206  columnNames.push_back( "x" ) ; columnTypes.push_back( "double" ) ;
207  columnNames.push_back( "y" ) ; columnTypes.push_back( "double" ) ;
208  columnNames.push_back( "z" ) ; columnTypes.push_back( "double" ) ;
209  columnNames.push_back( "sensitive" ) ; columnTypes.push_back( "bool" ) ;
210 
211  ITuple *vxdMap = mytuplefactory->create( mapName , mapTitle , columnNames , columnTypes ) ;
212 
213  // ladder and sensitive Detail
214  mapName = "vxdDetailMap" ;
215  mapTitle = "Detailed Map of part of the VTX" ;
216 
217  ITuple *vxdDetailMap = mytuplefactory->create( mapName , mapTitle , columnNames , columnTypes ) ;
218 
219  // map of points pointed on by VXDParametes::distanceToLadder()
220  mapName = "vxdPointsOnSurface" ;
221  mapTitle = "Map of Points On Surface" ;
222 
223  ITuple *vxdSurfaceMap = mytuplefactory->create( mapName , mapTitle , columnNames , columnTypes ) ;
224 
225  // map of point returned from intersection function between line and vxd
226  mapName = "vxdIntersectionPoints" ;
227  mapTitle = "Map of Intersection Points" ;
228 
229  ITuple *vxdIntersectionMap = mytuplefactory->create( mapName , mapTitle , columnNames , columnTypes ) ;
230 
231  // Histogram of Distances Surface to Ladder
232  IHistogram1D *distanceHisto = myhistofactory->createHistogram1D("distance",
233  100, 0, 0.00000000000001 ) ;
234 
235  // number of point to be tested
236  int nPoints = int( (testEnd-testStart)*(testEnd-testStart)*(zEnd-zStart)
237  / ( testStep * testStep * zStep )
238  + (testDetailEnd-testDetailStart)*(testDetailEnd-testDetailStart)*(zEnd-zStart)
239  / (testDetailStep * testDetailStep * zStep ) ) ;
240 
241  std::cout << "Testing " << nPoints << " Points to AIDA" << std::endl ;
242 
243  //counter for skipping
244  int skipped = 0 ;
245 
246  // first run over normal map
247  // x
248  for ( double x = testStart ; x < testEnd ; x += testStep ) {
249 
250  // Status
251  float curStatus = (x - testStart) / (testEnd - testStart) * 50 ;
252  if ( curStatus > ( myStatusP + 0.25 ) ) {
253  myStatusP = curStatus ;
254  std::cout << "." << std::flush ;
255  }
256  if ( curStatus > ( myStatusL + 10 ) ) {
257  myStatusL = curStatus ;
258  std::cout << floor(myStatusL) << "% finished" << std::endl ;
259  }
260 
261 
262  // y
263  for ( double y=testStart ; y < testEnd ; y += testStep ) {
264 
265  // z
266  for ( double z = zStart ; z < zEnd ; z += zStep ) {
267 
268  Vector3D p ( x , y , z ) ;
269  bool isPoint = vxdParams.isPointInLadder( p ) ;
270  bool isSensitive = vxdParams.isPointInSensitive( p ) ;
271 
272  if( isPoint ) nLadder ++;
273  if( isSensitive ) nSensitive ++;
274 
275  if( isPoint || isSensitive ) {
276 
277  // regular map
278  vxdMap->fill( 0 , p[0] ) ;
279  vxdMap->fill( 1 , p[1] ) ;
280  vxdMap->fill( 2 , p[2] ) ;
281  vxdMap->fill( 3 , isSensitive ) ;
282  vxdMap->addRow() ;
283 
284  }
285 
286  if( skipped >= skipDist ) {
287  // points pointed on
288  // ladder
289  Vector3D v = vxdParams.distanceToNearestLadder( p ) ;
290  Vector3D np ( p[0] + v[0] , p[1] + v[1] , p[2] + v[2] ) ;
291  // sensitive
292  v = vxdParams.distanceToNearestSensitive( p ) ;
293  Vector3D nps ( p[0] + v[0] , p[1] + v[1] , p[2] + v[2] ) ;
294 
295  vxdSurfaceMap->fill( 0 , np[0] ) ;
296  vxdSurfaceMap->fill( 1 , np[1] ) ;
297  vxdSurfaceMap->fill( 2 , np[2] ) ;
298  vxdSurfaceMap->fill( 3 , false ) ;
299  vxdSurfaceMap->addRow() ;
300 
301  vxdSurfaceMap->fill( 0 , nps[0] ) ;
302  vxdSurfaceMap->fill( 1 , nps[1] ) ;
303  vxdSurfaceMap->fill( 2 , nps[2] ) ;
304  vxdSurfaceMap->fill( 3 , true ) ;
305  vxdSurfaceMap->addRow() ;
306 
307  // distance histogram
308  v = vxdParams.distanceToNearestLadder( np ) ;
309  distanceHisto->fill( sqrt( v[0]*v[0] + v[1]*v[1] + v[2]*v[2] ) ) ;
310 
311  // reset counter
312  skipped = 0 ;
313  }
314  else {
315  skipped ++ ;
316  }
317 
318  } // z
319  } // y
320  } //x
321 
322 
323  // second run over detailed map
324 
325  // x
326  for ( double x = testDetailStart ; x < testDetailEnd ; x += testDetailStep ) {
327 
328  // Status
329  float curStatus = (x - testDetailStart) / (testDetailEnd - testDetailStart) * 50 + 50 ;
330  if ( curStatus > ( myStatusP + 0.25 ) ) {
331  myStatusP = curStatus ;
332  std::cout << "." << std::flush ;
333  }
334  if ( curStatus > ( myStatusL + 10 ) ) {
335  myStatusL = curStatus ;
336  std::cout << floor(myStatusL) << "% finished" << std::endl ;
337  }
338 
339 
340  // y
341  for ( double y=testDetailStart ; y < testDetailEnd ; y += testDetailStep ) {
342 
343  // z
344  for ( double z = zStart ; z < zEnd ; z += zStep ) {
345 
346  Vector3D p ( x , y , z ) ;
347  bool isPoint = vxdParams.isPointInLadder( p ) ;
348  bool isSensitive = vxdParams.isPointInSensitive( p ) ;
349 
350  if( isPoint ) nLadder ++;
351  if( isSensitive ) nSensitive ++;
352 
353  if( isPoint || isSensitive ) {
354 
355  // regular map
356  vxdDetailMap->fill( 0 , p[0] ) ;
357  vxdDetailMap->fill( 1 , p[1] ) ;
358  vxdDetailMap->fill( 2 , p[2] ) ;
359  vxdDetailMap->fill( 3 , isSensitive ) ;
360  vxdDetailMap->addRow() ;
361 
362  }
363 
364  } // z
365  } // y
366  } //x
367 
368 
369  // testing intersectionVXD
370 
371  std::cout << "Testing intersection " << std::flush ;
372 
373  Vector3D p0 (0. , 0. , 10. ) ;
374  Vector3D p1 (75. , 75. , 10. ) ;
375 
376  Vector3D lineVec ;
377 
378  for( double phi=-M_PI ; phi < M_PI ; phi += (M_PI/180) ) {
379  lineVec = Vector3D ( sin(phi) , cos(phi) , 0 ) ;
380  Vector3D iPoint0 = vxdParams.intersectionLadder( p0 , lineVec ) ;
381  Vector3D iPoint1 = vxdParams.intersectionLadder( p1 , lineVec ) ;
382  vxdIntersectionMap->fill( 0, iPoint0[0] ) ;
383  vxdIntersectionMap->fill( 1, iPoint0[1] ) ;
384  vxdIntersectionMap->fill( 2, iPoint0[2] ) ;
385  vxdIntersectionMap->fill( 3, false ) ;
386  vxdIntersectionMap->addRow() ;
387  vxdIntersectionMap->fill( 0, iPoint1[0] ) ;
388  vxdIntersectionMap->fill( 1, iPoint1[1] ) ;
389  vxdIntersectionMap->fill( 2, iPoint1[2] ) ;
390  vxdIntersectionMap->fill( 3, false ) ;
391  vxdIntersectionMap->addRow() ;
392  }
393 
394  mytree->commit() ;
395 
396  std::cout << "\ntotal tested points: " << nPoints << "\n"
397  << " in ladder: " << nLadder << " (" << nLadder/nPoints*100 << ") \n"
398  << " in sensitive: " << nSensitive << " (" << nSensitive/nPoints*100 << ")\n"
399  << "Points finished " << std::endl ;
400 
401 
402  // start visualizing
403 
404 
405 
406 }
407 
408 #endif
409 
410 
411 void testVXDPoint( const Vector3D p , const VXDParameters& vxdParams ) {
412 
413  std::cout << p[0] <<" " << p[1] <<" "<< p[2] ;
414  bool isPoint = vxdParams.isPointInLadder( p ) ;
415 
416  if( isPoint ) {
417  std::cout << p[0] <<" " << p[1] <<" "<< p[2] << " inside";
418  }
419 
420  std::cout << std::endl ;
421 
422 }
423 
424 void testVXDDist( const Vector3D p , const VXDParameters& vxdParams ) {
425 
426  Vector3D v = vxdParams.distanceToNearestLadder( p ) ;
427 
428  std::cout << "\nPoint (" << p[0] << "," << p[1] << "," << p[2] << ") "
429  << "\ndistance Vector " << v << std::endl;
430 // << "\ndistance Vector (" << v[0] << "," << v[1] << "," << v[2] << ") " << std::endl;
431 
432 
433  Vector3D paL = p + v ;
434 
435 // //FG: test the point found:
436 // std::cout << " is point in ladder : " << vxdParams.isPointInLadder( paL ) << std::endl ;
437 
438  // how far are we still from the ladder ?
439  Vector3D v2 = vxdParams.distanceToNearestLadder( paL ) ;
440 
441  std::cout << " distance to ladder : " << v2.r() << std::endl ;
442 
443 
444 
445 
446 // bool isFound = false ;
447 
448 // double myStep =1 ;
449 // double lBound = 1 ;
450 // double uBound = 1 ;
451 
452 // for(double r1 = lBound ; r1 <= uBound ; r1+= myStep ) {
453 // for(double r2 = lBound ; r2 <= uBound ; r2+= myStep ) {
454 // for(double r3 = lBound ; r3 <= uBound ; r3+= myStep ) {
455 
456 // Vector3D iP ( p[0] + r1* v[0] , p[1] + r2* v[1] , p[2] + r3*v[2] ) ;
457 
458 // bool isPoint = vxdParams.isPointInLadder( iP ) ;
459 // if( isPoint ) {
460 // std::cout << "Vector point withhin Ladder at r =" << r1 <<" , " << r2 <<" , " << r3 << std::endl ;
461 // std::cout << " " << iP[0] << " , " << iP[1] << " , " << iP[2] << std::endl ;
462 // isFound = true ;
463 // break ;
464 // }
465 // else{
466 // std::cout << "ERROR Vector NOT pointin in Ladder." << std::endl ;
467 // }
468 // }
469 // if (isFound ) break ;
470 // }
471 // if (isFound ) break ;
472 // }
473 
474 
475 }
476 
477 void testAllVXDPoint( const VXDParameters& vxdParams ) {
478 
479  double testStart = -65 ;
480  double testEnd = 65 ;
481  double testStep = 0.05 ;
482  double zStart = 10 ;
483  double zEnd = 11 ;
484  double zStep = 1 ;
485 
486  float myStatusP = 0, myStatusL = 0;
487  int nLadder = 0 , nSensitive = 0 ;
488 
489  std::ofstream myfile;
490  const char *fileName = "./gearVXDTest.txt" ;
491  myfile.open (fileName, std::ios::out | std::ios::trunc );
492 
493 
494 #ifdef GEAR_USE_AIDA
495 // IHistogram2D *ladderMap = myhistofactory->createHistogram2D("ladderMap",
496 // 500, testStart, testEnd,
497 // 500, testStart, testEnd );
498  std::string mapName = "vxdMap" ;
499  std::string mapTitle = "Map of the VTX" ;
500  std::vector<std::string> columnNames ;
501  std::vector<std::string> columnTypes ;
502 
503  columnNames.push_back( "x" ) ; columnTypes.push_back( "double" ) ;
504  columnNames.push_back( "y" ) ; columnTypes.push_back( "double" ) ;
505  columnNames.push_back( "z" ) ; columnTypes.push_back( "double" ) ;
506  columnNames.push_back( "sensitive" ) ; columnTypes.push_back( "bool" ) ;
507 
508  ITuple *ladderMap = mytuplefactory->create( mapName , mapTitle , columnNames , columnTypes ) ;
509 
510 
511 #endif
512 
513  int nPoints = int( (testEnd-testStart)*(testEnd-testStart)*(zEnd-zStart)
514  / ( testStep * testStep * zStep ) ) ;
515  std::cout << "Testing " << nPoints << " Points to AIDA" << std::endl ;
516 
517  // x
518  for ( double x=testStart ; x < testEnd ; x += testStep ) {
519 
520  // Status
521  float curStatus = (x - testStart) / (testEnd - testStart) * 100 ;
522  if ( curStatus > ( myStatusP + 0.25 ) ) {
523  myStatusP = curStatus ;
524  std::cout << "." << std::flush ;
525  }
526  if ( curStatus > ( myStatusL + 10 ) ) {
527  myStatusL = curStatus ;
528  std::cout << floor(myStatusL) << "% finished" << std::endl ;
529  }
530 
531  // y
532  for ( double y=testStart ; y < testEnd ; y += testStep ) {
533 
534  // z
535  for ( double z = zStart ; z < zEnd ; z += zStep ) {
536 
537  Vector3D p ( x , y , z ) ;
538  bool isPoint = vxdParams.isPointInLadder( p ) ;
539  bool isSensitive = vxdParams.isPointInSensitive( p ) ;
540 
541  if( isPoint ) nLadder ++;
542  if( isSensitive ) nSensitive ++;
543 
544  if( isPoint || isSensitive ) {
545  myfile << p[0] <<" " << p[1] <<" "<< p[2] << " " << isSensitive << "\n" ;
546 
547 #ifdef GEAR_USE_AIDA
548  ladderMap->fill( 0 , p[0] ) ;
549  ladderMap->fill( 1 , p[1] ) ;
550  ladderMap->fill( 2 , p[2] ) ;
551  ladderMap->fill( 3 , isSensitive ) ;
552  ladderMap->addRow() ;
553 #endif
554  }
555 
556  } // z
557  } // y
558  } //x
559 
560  myfile.close() ;
561 
562  std::cout << "\ntotal tested points: " << nPoints << "\n"
563  << " in ladder: " << nLadder << " (" << nLadder/nPoints*100 << ") \n"
564  << " in sensitive: " << nSensitive << " (" << nSensitive/nPoints*100 << ")\n"
565  << "Points finished " << std::endl ;
566 
567 }
568 
569 
570 void testAllVXDDist( const VXDParameters& vxdParams ) {
571 
572  double testStart = -65 ;
573  double testEnd = 65 ;
574  double testStep = 1 ;
575  double zStart = -120 ;
576  double zEnd = 120 ;
577  double zStep = 10 ;
578 
579  float myStatusP = 0, myStatusL = 0;
580 
581  std::ofstream myfile, myHisto;
582  const char *fileName = "./gearVXDTest_out.txt" ;
583  const char *histoName = "./gearVXDTest_Histogram.txt" ;
584  myfile.open (fileName, std::ios::out | std::ios::trunc );
585  myHisto.open( histoName, std::ios::out | std::ios::trunc ) ;
586 
587  int nPoints = int ( (testEnd-testStart)*(testEnd-testStart)*(zEnd-zStart)
588  / ( testStep * testStep * zStep ) ) ;
589  std::cout << "Testing " << nPoints << " Points to file : " << fileName <<std::endl ;
590 
591  // x
592  for ( double x=testStart ; x < testEnd ; x += testStep ) {
593 
594  // Status
595  float curStatus = (x - testStart) / (testEnd - testStart) * 100 ;
596  if ( curStatus > ( myStatusP + 0.25 ) ) {
597  myStatusP = curStatus ;
598  std::cout << "." << std::flush ;
599  }
600  if ( curStatus > ( myStatusL + 10 ) ) {
601  myStatusL = curStatus ;
602  std::cout << floor(myStatusL) << "% finished" << std::endl ;
603  }
604 
605  // y
606  for ( double y=testStart ; y < testEnd ; y += testStep ) {
607 
608  // z
609  for ( double z = zStart ; z < zEnd ; z += zStep ) {
610 
611  Vector3D p ( x , y , z ) ;
612  Vector3D nv = vxdParams.distanceToNearestLadder( p ) ;
613  Vector3D np ( p[0] + nv[0] , p[1] + nv[1] , p[2] + nv[2] ) ;
614 
615  //bool isPoint = vxdParams.isPointInLadder( np ) ;
616 
617  //if( !isPoint ) {
618  Vector3D nv2 = vxdParams.distanceToNearestLadder( np ) ;
619  myHisto << sqrt( nv2[0]*nv2[0] + nv2[1]*nv2[1] + nv2[2]*nv2[2] ) << "\n" ;
620  myfile << np[0] <<" " << np[1] <<" "<< np[2] << 0 << "\n" ;
621  //}
622 
623  } // z
624  } // y
625  } //x
626 
627  myfile.close() ;
628  myHisto.close() ;
629 
630  std::cout << "Points finished " << std::endl ;
631 }
virtual bool isPointInSensitive(Vector3D p, SensorID *sensorID=0) const =0
returns wheter a point is inside a sensitive volume- if sensorID != 0 the sensorID is returned in cas...
Geometry properties of a vertex detector needed for reconstruction code.
virtual const ZPlanarParameters & getVXDParameters() const =0
Get the VXD parameters.
virtual bool isPointInLadder(Vector3D p) const =0
returns whether a point is inside a ladder
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 Vector3D distanceToNearestLadder(Vector3D p) const =0
returns vector from point to nearest ladder
virtual Vector3D intersectionLadder(Vector3D p, Vector3D v) const =0
returns the first point where a given straight line (parameters point p and direction v) crosses a la...
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 Vector3D distanceToNearestSensitive(Vector3D p) const =0
returns vector from point to nearest sensitive volume