MagneticField/interface/Field.h

0001 #ifndef SimG4Core_Field_H
0002 #define SimG4Core_Field_H
0003
0004 #include "G4MagneticField.hh"
0005
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007
0008 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0009
0010 #include <CLHEP/Units/SystemOfUnits.h>
0011
0012 namespace sim {
0013   class Field final : public G4MagneticField {
0014   public:
0015     Field(const MagneticField *f, double d);
0016     ~Field() override;
0017     inline void GetFieldValue(const G4double p[], G4double b[3]) const override;
0018
0019   private:
0020     const MagneticField *theCMSMagneticField;
0021     double theDelta;
0022
0023     mutable double oldx[3];
0024     mutable double oldb[3];
0025   };
0026 };  // namespace sim
0027
0028 void sim::Field::GetFieldValue(const G4double xyz[], G4double bfield[3]) const {
0029   if (std::abs(oldx[0] - xyz[0]) > theDelta || std::abs(oldx[1] - xyz[1]) > theDelta ||
0030       std::abs(oldx[2] - xyz[2]) > theDelta) {
0031     constexpr float lunit = (1.0 / CLHEP::cm);
0032     GlobalPoint ggg((float)(xyz[0]) * lunit, (float)(xyz[1]) * lunit, (float)(xyz[2]) * lunit);
0033     GlobalVector v = theCMSMagneticField->inTesla(ggg);
0034
0035     constexpr float btesla = CLHEP::tesla;
0036     oldb[0] = (v.x() * btesla);
0037     oldb[1] = (v.y() * btesla);
0038     oldb[2] = (v.z() * btesla);
0039     oldx[0] = xyz[0];
0040     oldx[1] = xyz[1];
0041     oldx[2] = xyz[2];
0042   }
0043
0044   bfield[0] = oldb[0];
0045   bfield[1] = oldb[1];
0046   bfield[2] = oldb[2];
0047 }
0048
0049 #endif