ParticlePropagator/src/MagneticFieldMap.cc

0001 #include "MagneticField/Engine/interface/MagneticField.h"
0002 #include "FastSimulation/ParticlePropagator/interface/MagneticFieldMap.h"
0003 #include "FastSimulation/TrackerSetup/interface/TrackerInteractionGeometry.h"
0004
0005 #include <iostream>
0006
0007 MagneticFieldMap::MagneticFieldMap(const MagneticField* pMF, const TrackerInteractionGeometry* myGeo)
0008     : pMF_(pMF),
0009       geometry_(myGeo),
0010       bins(101),
0011       fieldBarrelHistos(200, static_cast<std::vector<double> >(std::vector<double>(bins, static_cast<double>(0.)))),
0012       fieldEndcapHistos(200, static_cast<std::vector<double> >(std::vector<double>(bins, static_cast<double>(0.)))),
0013       fieldBarrelBinWidth(200, static_cast<double>(0.)),
0014       fieldBarrelZMin(200, static_cast<double>(0.)),
0015       fieldEndcapBinWidth(200, static_cast<double>(0.)),
0016       fieldEndcapRMin(200, static_cast<double>(0.))
0017
0018 {
0019   std::list<TrackerLayer>::const_iterator cyliter;
0020   std::list<TrackerLayer>::const_iterator cylitBeg = geometry_->cylinderBegin();
0021   std::list<TrackerLayer>::const_iterator cylitEnd = geometry_->cylinderEnd();
0022
0023   // Prepare the histograms
0024   // std::cout << "Prepare magnetic field local database for FAMOS speed-up" << std::endl;
0025   for (cyliter = cylitBeg; cyliter != cylitEnd; ++cyliter) {
0026     int layer = cyliter->layerNumber();
0027     //    cout << " Fill Histogram " << hist << endl;
0028
0029     // Cylinder bounds
0030     double zmin = 0.;
0031     double zmax;
0032     double rmin = 0.;
0033     double rmax;
0034     if (cyliter->forward()) {
0035       zmax = cyliter->disk()->position().z();
0036       rmax = cyliter->disk()->outerRadius();
0037     } else {
0038       zmax = cyliter->cylinder()->bounds().length() / 2.;
0039       rmax = cyliter->cylinder()->bounds().width() / 2. - cyliter->cylinder()->bounds().thickness() / 2.;
0040     }
0041
0042     // Histograms
0043     double step;
0044
0045     // Disk histogram characteristics
0046     step = (rmax - rmin) / (bins - 1);
0047     fieldEndcapBinWidth[layer] = step;
0048     fieldEndcapRMin[layer] = rmin;
0049
0050     // Fill the histo
0051     int endcapBin = 0;
0052     for (double radius = rmin + step / 2.; radius < rmax + step; radius += step) {
0053       double field = inTeslaZ(GlobalPoint(radius, 0., zmax));
0054       fieldEndcapHistos[layer][endcapBin++] = field;
0055     }
0056
0057     // Barrel Histogram characteritics
0058     step = (zmax - zmin) / (bins - 1);
0059     fieldBarrelBinWidth[layer] = step;
0060     fieldBarrelZMin[layer] = zmin;
0061
0062     // Fill the histo
0063     int barrelBin = 0;
0064     for (double zed = zmin + step / 2.; zed < zmax + step; zed += step) {
0065       double field = inTeslaZ(GlobalPoint(rmax, 0., zed));
0066       fieldBarrelHistos[layer][barrelBin++] = field;
0067     }
0068   }
0069 }
0070
0071 const GlobalVector MagneticFieldMap::inTesla(const GlobalPoint& gp) const {
0072   if (!pMF_) {
0073     return GlobalVector(0., 0., 4.);
0074   } else {
0075     return pMF_->inTesla(gp);
0076   }
0077 }
0078
0079 const GlobalVector MagneticFieldMap::inTesla(const TrackerLayer& aLayer, double coord, int success) const {
0080   if (!pMF_) {
0081     return GlobalVector(0., 0., 4.);
0082   } else {
0083     return GlobalVector(0., 0., inTeslaZ(aLayer, coord, success));
0084   }
0085 }
0086
0087 const GlobalVector MagneticFieldMap::inKGauss(const GlobalPoint& gp) const { return inTesla(gp) * 10.; }
0088
0089 const GlobalVector MagneticFieldMap::inInverseGeV(const GlobalPoint& gp) const { return inKGauss(gp) * 2.99792458e-4; }
0090
0091 double MagneticFieldMap::inTeslaZ(const GlobalPoint& gp) const { return pMF_ ? pMF_->inTesla(gp).z() : 4.0; }
0092
0093 double MagneticFieldMap::inTeslaZ(const TrackerLayer& aLayer, double coord, int success) const {
0094   if (!pMF_) {
0095     return 4.;
0096   } else {
0097     // Find the relevant histo
0098     double theBinWidth;
0099     double theXMin;
0100     unsigned layer = aLayer.layerNumber();
0101     const std::vector<double>* theHisto;
0102
0103     if (success == 1) {
0104       theHisto = theFieldBarrelHisto(layer);
0105       theBinWidth = fieldBarrelBinWidth[layer];
0106       theXMin = fieldBarrelZMin[layer];
0107     } else {
0108       theHisto = theFieldEndcapHisto(layer);
0109       theBinWidth = fieldEndcapBinWidth[layer];
0110       theXMin = fieldEndcapRMin[layer];
0111     }
0112
0113     // Find the relevant bin
0114     double x = fabs(coord);
0115     unsigned bin = (unsigned)((x - theXMin) / theBinWidth);
0116     if (bin + 1 == (unsigned)bins)
0117       bin -= 1;  // Add a protection against coordinates near the layer edge
0118     double x1 = theXMin + (bin - 0.5) * theBinWidth;
0119     double x2 = x1 + theBinWidth;
0120
0121     // Determine the field
0122     double field1 = (*theHisto)[bin];
0123     double field2 = (*theHisto)[bin + 1];
0124
0125     return field1 + (field2 - field1) * (x - x1) / (x2 - x1);
0126   }
0127 }
0128
0129 double MagneticFieldMap::inKGaussZ(const GlobalPoint& gp) const { return inTeslaZ(gp) / 10.; }
0130
0131 double MagneticFieldMap::inInverseGeVZ(const GlobalPoint& gp) const { return inKGaussZ(gp) * 2.99792458e-4; }