/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/FastSimulation/ParticlePropagator/src/MagneticFieldMap.cc

Go to the documentation of this file.

#include "MagneticField/Engine/interface/MagneticField.h"
#include "FastSimulation/ParticlePropagator/interface/MagneticFieldMap.h"
#include "FastSimulation/TrackerSetup/interface/TrackerInteractionGeometry.h"

#include <iostream>

MagneticFieldMap::MagneticFieldMap(const MagneticField* pMF,
                                   const TrackerInteractionGeometry* myGeo) : 
  pMF_(pMF), 
  geometry_(myGeo),
  bins(101), 
  fieldBarrelHistos(200,static_cast<std::vector<double> >
                    (std::vector<double>(bins,static_cast<double>(0.)))),
  fieldEndcapHistos(200,static_cast<std::vector<double> > 
                     (std::vector<double>(bins,static_cast<double>(0.)))),
  fieldBarrelBinWidth(200,static_cast<double>(0.)),
  fieldBarrelZMin(200,static_cast<double>(0.)),
  fieldEndcapBinWidth(200,static_cast<double>(0.)),
  fieldEndcapRMin(200,static_cast<double>(0.))

{
  
  std::list<TrackerLayer>::const_iterator cyliter;
  std::list<TrackerLayer>::const_iterator cylitBeg=geometry_->cylinderBegin();
  std::list<TrackerLayer>::const_iterator cylitEnd=geometry_->cylinderEnd();
  
  // Prepare the histograms
  // std::cout << "Prepare magnetic field local database for FAMOS speed-up" << std::endl;
  for ( cyliter=cylitBeg; cyliter != cylitEnd; ++cyliter ) {
    int layer = cyliter->layerNumber();
    //    cout << " Fill Histogram " << hist << endl;

    // Cylinder bounds
    double zmin = 0.;
    double zmax; 
    double rmin = 0.;
    double rmax; 
    if ( cyliter->forward() ) {
      zmax = cyliter->disk()->position().z();
      rmax = cyliter->disk()->outerRadius();
    } else {
      zmax = cyliter->cylinder()->bounds().length()/2.;
      rmax = cyliter->cylinder()->bounds().width()/2.
           - cyliter->cylinder()->bounds().thickness()/2.;
    }

    // Histograms
    double step;

    // Disk histogram characteristics
    step = (rmax-rmin)/(bins-1);
    fieldEndcapBinWidth[layer] = step;
    fieldEndcapRMin[layer] = rmin;

    // Fill the histo
    int endcapBin = 0;
    for ( double radius=rmin+step/2.; radius<rmax+step; radius+=step ) {
      double field = inTeslaZ(GlobalPoint(radius,0.,zmax));
      fieldEndcapHistos[layer][endcapBin++] = field;
    }

    // Barrel Histogram characteritics
    step = (zmax-zmin)/(bins-1);
    fieldBarrelBinWidth[layer] = step;
    fieldBarrelZMin[layer] = zmin;

    // Fill the histo
    int barrelBin = 0;
    for ( double zed=zmin+step/2.; zed<zmax+step; zed+=step ) {
      double field = inTeslaZ(GlobalPoint(rmax,0.,zed));
      fieldBarrelHistos[layer][barrelBin++] = field;
    }
  }
}


const GlobalVector
MagneticFieldMap::inTesla( const GlobalPoint& gp) const {

  if (!pMF_) {
    return GlobalVector( 0., 0., 4.);
  } else {
    return pMF_->inTesla(gp);
  }

}

const GlobalVector
MagneticFieldMap::inTesla(const TrackerLayer& aLayer, double coord, int success) const {

  if (!pMF_) {
    return GlobalVector( 0., 0., 4.);
  } else {
    return GlobalVector(0.,0.,inTeslaZ(aLayer,coord,success));
  }

}

const GlobalVector 
MagneticFieldMap::inKGauss( const GlobalPoint& gp) const {
  
  return inTesla(gp) * 10.;

}

const GlobalVector
MagneticFieldMap::inInverseGeV( const GlobalPoint& gp) const {
  
  return inKGauss(gp) * 2.99792458e-4;

} 

double 
MagneticFieldMap::inTeslaZ(const GlobalPoint& gp) const {

    return pMF_ ? pMF_->inTesla(gp).z() : 4.0;

}

double
MagneticFieldMap::inTeslaZ(const TrackerLayer& aLayer, double coord, int success) const 
{

  if (!pMF_) {
    return 4.;
  } else {
    // Find the relevant histo
    double theBinWidth;
    double theXMin;
    unsigned layer = aLayer.layerNumber();
    const std::vector<double>* theHisto;

    if ( success == 1 ) { 
      theHisto = theFieldBarrelHisto(layer);
      theBinWidth = fieldBarrelBinWidth[layer];
      theXMin = fieldBarrelZMin[layer];
    } else {
      theHisto = theFieldEndcapHisto(layer);
      theBinWidth = fieldEndcapBinWidth[layer];
      theXMin = fieldEndcapRMin[layer];
    }
    
    // Find the relevant bin
    double x = fabs(coord);
    unsigned bin = (unsigned) ((x-theXMin)/theBinWidth);
    if ( bin+1 == (unsigned)bins ) bin -= 1; // Add a protection against coordinates near the layer edge
    double x1 = theXMin + (bin-0.5)*theBinWidth;
    double x2 = x1+theBinWidth;      

    // Determine the field
    double field1 = (*theHisto)[bin];
    double field2 = (*theHisto)[bin+1];

    return field1 + (field2-field1) * (x-x1)/(x2-x1);

  }

}

double 
MagneticFieldMap::inKGaussZ(const GlobalPoint& gp) const {

    return inTeslaZ(gp)/10.;

}

double 
MagneticFieldMap::inInverseGeVZ(const GlobalPoint& gp) const {

   return inKGaussZ(gp) * 2.99792458e-4;

}