MultipleScatteringParametrisation.cc

Go to the documentation of this file.

#include <vector>
#include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisation.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "RecoTracker/TkMSParametrization/interface/PixelRecoPointRZ.h"
#include "RecoTracker/TkMSParametrization/interface/PixelRecoLineRZ.h"

template <class T>
inline T sqr(T t) {
  return t * t;
}

#include "MSLayersKeeper.h"
#include "MSLayersKeeperX0AtEta.h"
#include "MSLayersKeeperX0Averaged.h"
#include "MSLayersKeeperX0DetLayer.h"
#include "MSLayersAtAngle.h"

//#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"

#include <iostream>

using namespace std;

const float MultipleScatteringParametrisation::x0ToSigma = 0.0136f;

namespace {
  struct Keepers {
    MSLayersKeeperX0DetLayer x0DetLayer;
    MSLayersKeeperX0AtEta x0AtEta;
    MSLayersKeeperX0Averaged x0Averaged;
    MSLayersKeeper *keepers[3];  // {&x0DetLayer,&x0AtEta,&x0Averaged};
    bool isInitialised;          // =false;
    MSLayersKeeper const *operator()(int i) const { return keepers[i]; }
    void init(const edm::EventSetup &iSetup) {
      if (isInitialised)
        return;
      for (auto x : keepers)
        x->init(iSetup);
      isInitialised = true;
    }
    Keepers() : keepers{&x0DetLayer, &x0AtEta, &x0Averaged}, isInitialised(false) {}
  };

  thread_local Keepers keepers;
  //NOTE: This is being used to globally cache information from the EventSetup
  // this should not be done so we need this code changed.
  //NOTE; the thread_local only works in this case because MultipleScateringParametrisation
  // instances are only ever created on the stack and not the heap.

}  // namespace

void MultipleScatteringParametrisation::initKeepers(const edm::EventSetup &iSetup) { keepers.init(iSetup); }

//using namespace PixelRecoUtilities;
//----------------------------------------------------------------------
void MultipleScatteringParametrisation::init(const DetLayer *layer, const edm::EventSetup &iSetup, X0Source x0Source) {
  theLayerKeeper = keepers(x0Source);

  // FIXME not thread safe: move elsewhere...
  initKeepers(iSetup);

  if (!layer)
    return;
  theLayer = theLayerKeeper->layer(layer);
}

//----------------------------------------------------------------------
float MultipleScatteringParametrisation::operator()(float pT, float cotTheta, float) const {
  float sumX0D = theLayer.sumX0D(cotTheta);
  return x0ToSigma * sumX0D / pT;
}

//----------------------------------------------------------------------
float MultipleScatteringParametrisation::operator()(float pT,
                                                    float cotTheta,
                                                    const PixelRecoPointRZ &pointI,
                                                    float tip) const {
  PixelRecoLineRZ lineIO(pointI, cotTheta, tip);
  PixelRecoPointRZ pointO = theLayer.crossing(lineIO).first;

  const MSLayersAtAngle &layersAtEta = theLayerKeeper->layers(cotTheta);

  float sumX0D = layersAtEta.sumX0D(pointI, pointO);
  return x0ToSigma * sumX0D / pT;
}

float MultipleScatteringParametrisation::operator()(float pT,
                                                    float cotTheta,
                                                    const PixelRecoPointRZ &pointI,
                                                    int il) const {
  PixelRecoLineRZ lineIO(pointI, cotTheta);
  PixelRecoPointRZ pointO = theLayer.crossing(lineIO).first;

  const MSLayersAtAngle &layersAtEta = theLayerKeeper->layers(cotTheta);

  float sumX0D = layersAtEta.sumX0D(il, theLayer.seqNum(), pointI, pointO);
  return x0ToSigma * sumX0D / pT;
}

//----------------------------------------------------------------------
float MultipleScatteringParametrisation::operator()(float pT,
                                                    const PixelRecoPointRZ &pointI,
                                                    const PixelRecoPointRZ &pointO,
                                                    Consecutive consecutive,
                                                    float tip) const {
  PixelRecoLineRZ lineIO(pointI, pointO, tip);
  PixelRecoPointRZ pointM = theLayer.crossing(lineIO).first;
  float cotTheta = lineIO.cotLine();

  if (consecutive == useConsecutive) {
    float dist = fabs((pointO.r() - pointM.r()) * (pointM.r() - pointI.r()) / (pointO.r() - pointI.r()));
    return x0ToSigma * sqrt(theLayer.x0(cotTheta)) * dist / pT;
  } else {
    const MSLayersAtAngle &layersAtEta = theLayerKeeper->layers(cotTheta);
    float sumX0D = layersAtEta.sumX0D(pointI, pointM, pointO);
    return x0ToSigma * sumX0D / pT;
  }
}

float MultipleScatteringParametrisation::operator()(float pT,
                                                    const PixelRecoPointRZ &pointV,
                                                    const PixelRecoPointRZ &pointO,
                                                    int ol) const {
  PixelRecoLineRZ lineIO(pointV, pointO);
  PixelRecoPointRZ pointI = theLayer.crossing(lineIO).first;
  float cotTheta = lineIO.cotLine();

  const MSLayersAtAngle &layersAtEta = theLayerKeeper->layers(cotTheta);
  float sumX0D = layersAtEta.sumX0D(pointV.z(), theLayer.seqNum(), ol, pointI, pointO);
  return x0ToSigma * sumX0D / pT;
}