File indexing completed on 2023-03-17 11:26:26
0001
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 }
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 }