GsfTracking/test/GsfMaterialEffectsUpdator_t.cpp

0001 // test of  GsfMaterialEffectsUpdator
0002
0003 #include "TrackingTools/GsfTracking/interface/GsfMaterialEffectsUpdator.h"
0004 #include "TrackingTools/MaterialEffects/interface/MultipleScatteringUpdator.h"
0005 #include "TrackingTools/GsfTracking/interface/GsfCombinedMaterialEffectsUpdator.h"
0006 #include "TrackingTools/GsfTracking/interface/GsfMaterialEffectsAdapter.h"
0007
0008 #include "TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h"
0009
0010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0011 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0012 #include "DataFormats/GeometrySurface/interface/Surface.h"
0013
0014 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0015
0016 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0017 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0018 #include "MagneticField/Engine/interface/MagneticField.h"
0019
0020 namespace {
0021
0022   struct M5T : public MagneticField {
0023     M5T() : m(0., 0., 5.) {}
0024     virtual GlobalVector inTesla(const GlobalPoint&) const { return m; }
0025
0026     GlobalVector m;
0027   };
0028
0029 }  // namespace
0030
0031 #include "FWCore/Utilities/interface/HRRealTime.h"
0032 #include <iostream>
0033 #include <vector>
0034
0035 void st() {}
0036 void en() {}
0037
0038 int main(int argc, char* arg[]) {
0039   std::string file("BetheHeitler_cdfmom_nC6_O5.par");
0040   if (argc < 2) {
0041     std::cerr << "parameter file not given in input: default used" << std::endl;
0042   } else
0043     file = arg[1];
0044
0045   std::cout << "using file " << file << std::endl;
0046
0047   GsfBetheHeitlerUpdator bhu(file, 0);
0048   GsfBetheHeitlerUpdator bhu1(file, 1);
0049   GsfBetheHeitlerUpdator bhu2(file, 2);
0050   GsfMaterialEffectsAdapter msu(MultipleScatteringUpdator(bhu.mass()));
0051
0052   GsfCombinedMaterialEffectsUpdator comb(msu, bhu);
0053
0054   GsfMaterialEffectsUpdator* meus[] = {&msu, &bhu, &bhu1, &bhu2, &comb};
0055
0056   double neverKnow = 0;
0057   for (int j = 0; j != 5; ++j) {
0058     GsfMaterialEffectsUpdator* meu = meus[j];
0059
0060     Basic3DVector<float> axis(0.5, 1., 1);
0061
0062     Surface::RotationType rot(axis, 0.5 * M_PI);
0063
0064     Surface::PositionType pos(0., 0., 0.);
0065
0066     Plane::PlanePointer plane = Plane::build(pos, rot, new RectangularPlaneBounds(1., 1., 1));
0067     plane->setMediumProperties(MediumProperties(0.1, 0.3));
0068     M5T const m;
0069
0070     edm::HRTimeDiffType totT = 0;
0071     int n = 0;
0072     bool printIt = true;
0073     for (int i = 0; i != 100000; ++i) {
0074       LocalTrajectoryParameters tp(1. / (10. + 0.01 * i), 1., 1., 0., 0., 1.);
0075       LocalTrajectoryError lerr(1., 1., 0.1, 0.1, 0.1);
0076
0077       TrajectoryStateOnSurface tsos(tp, lerr, *plane, &m, SurfaceSideDefinition::beforeSurface);
0078       if (printIt) {
0079         std::cout << tsos.globalMomentum() << std::endl;
0080         std::cout << tsos.localError().matrix() << std::endl;
0081         std::cout << tsos.weight() << std::endl;
0082       }
0083       st();
0084       totT -= edm::hrRealTime();
0085       TrajectoryStateOnSurface tsos2 = meu->updateState(tsos, alongMomentum);
0086       totT += edm::hrRealTime();
0087       ++n;
0088       en();
0089       neverKnow += tsos2.globalMomentum().perp();
0090       if (printIt) {
0091         std::cout << tsos2.globalMomentum() << std::endl;
0092         std::cout << tsos2.localError().matrix() << std::endl;
0093         std::cout << tsos2.weight() << std::endl;
0094         printIt = false;
0095       }
0096     }
0097
0098     std::cout << "\nupdate  time " << double(totT) / double(n) << std::endl;
0099   }
0100   return neverKnow != 0 ? 0 : 20;
0101 }