ShallowGainCalibration.cc

Go to the documentation of this file.

#include "CalibTracker/SiStripCommon/interface/ShallowGainCalibration.h"
#include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
#include "CalibTracker/Records/interface/SiStripGainRcd.h"

using namespace edm;
using namespace reco;
using namespace std;

ShallowGainCalibration::ShallowGainCalibration(const edm::ParameterSet& iConfig)
    : tracks_token_(consumes<edm::View<reco::Track>>(iConfig.getParameter<edm::InputTag>("Tracks"))),
      association_token_(consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("Tracks"))),
      Suffix(iConfig.getParameter<std::string>("Suffix")),
      Prefix(iConfig.getParameter<std::string>("Prefix")) {
  produces<std::vector<int>>(Prefix + "trackindex" + Suffix);
  produces<std::vector<unsigned int>>(Prefix + "rawid" + Suffix);
  produces<std::vector<double>>(Prefix + "localdirx" + Suffix);
  produces<std::vector<double>>(Prefix + "localdiry" + Suffix);
  produces<std::vector<double>>(Prefix + "localdirz" + Suffix);
  produces<std::vector<unsigned short>>(Prefix + "firststrip" + Suffix);
  produces<std::vector<unsigned short>>(Prefix + "nstrips" + Suffix);
  produces<std::vector<bool>>(Prefix + "saturation" + Suffix);
  produces<std::vector<bool>>(Prefix + "overlapping" + Suffix);
  produces<std::vector<bool>>(Prefix + "farfromedge" + Suffix);
  produces<std::vector<unsigned int>>(Prefix + "charge" + Suffix);
  produces<std::vector<double>>(Prefix + "path" + Suffix);
#ifdef ExtendedCALIBTree
  produces<std::vector<double>>(Prefix + "chargeoverpath" + Suffix);
#endif
  produces<std::vector<unsigned char>>(Prefix + "amplitude" + Suffix);
  produces<std::vector<double>>(Prefix + "gainused" + Suffix);
  produces<std::vector<double>>(Prefix + "gainusedTick" + Suffix);
}

void ShallowGainCalibration::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  auto trackindex = std::make_unique<std::vector<int>>();
  auto rawid = std::make_unique<std::vector<unsigned int>>();
  auto localdirx = std::make_unique<std::vector<double>>();
  auto localdiry = std::make_unique<std::vector<double>>();
  auto localdirz = std::make_unique<std::vector<double>>();
  auto firststrip = std::make_unique<std::vector<unsigned short>>();
  auto nstrips = std::make_unique<std::vector<unsigned short>>();
  auto saturation = std::make_unique<std::vector<bool>>();
  auto overlapping = std::make_unique<std::vector<bool>>();
  auto farfromedge = std::make_unique<std::vector<bool>>();
  auto charge = std::make_unique<std::vector<unsigned int>>();
  auto path = std::make_unique<std::vector<double>>();
#ifdef ExtendedCALIBTree
  auto chargeoverpath = std::make_unique<std::vector<double>>();
#endif
  auto amplitude = std::make_unique<std::vector<unsigned char>>();
  auto gainused = std::make_unique<std::vector<double>>();
  auto gainusedTick = std::make_unique<std::vector<double>>();

  edm::ESHandle<TrackerGeometry> theTrackerGeometry;
  iSetup.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
  m_tracker = &(*theTrackerGeometry);
  edm::ESHandle<SiStripGain> gainHandle;
  iSetup.get<SiStripGainRcd>().get(gainHandle);
  edm::Handle<edm::View<reco::Track>> tracks;
  iEvent.getByToken(tracks_token_, tracks);
  edm::Handle<TrajTrackAssociationCollection> associations;
  iEvent.getByToken(association_token_, associations);

  for (TrajTrackAssociationCollection::const_iterator association = associations->begin();
       association != associations->end();
       association++) {
    const Trajectory* traj = association->key.get();
    const reco::Track* track = association->val.get();

    vector<TrajectoryMeasurement> measurements = traj->measurements();
    for (vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin();
         measurement_it != measurements.end();
         measurement_it++) {
      TrajectoryStateOnSurface trajState = measurement_it->updatedState();
      if (!trajState.isValid())
        continue;

      const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
      const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
      const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
      const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
      const SiPixelRecHit* sipixelhit = dynamic_cast<const SiPixelRecHit*>(hit);

      const SiPixelCluster* PixelCluster = nullptr;
      const SiStripCluster* StripCluster = nullptr;
      uint32_t DetId = 0;

      for (unsigned int h = 0; h < 2; h++) {
        if (!sistripmatchedhit && h == 1) {
          continue;
        } else if (sistripmatchedhit && h == 0) {
          StripCluster = &sistripmatchedhit->monoCluster();
          DetId = sistripmatchedhit->monoId();
        } else if (sistripmatchedhit && h == 1) {
          StripCluster = &sistripmatchedhit->stereoCluster();
          ;
          DetId = sistripmatchedhit->stereoId();
        } else if (sistripsimplehit) {
          StripCluster = (sistripsimplehit->cluster()).get();
          DetId = sistripsimplehit->geographicalId().rawId();
        } else if (sistripsimple1dhit) {
          StripCluster = (sistripsimple1dhit->cluster()).get();
          DetId = sistripsimple1dhit->geographicalId().rawId();
        } else if (sipixelhit) {
          PixelCluster = (sipixelhit->cluster()).get();
          DetId = sipixelhit->geographicalId().rawId();
        } else {
          continue;
        }

        LocalVector trackDirection = trajState.localDirection();
        double cosine = trackDirection.z() / trackDirection.mag();
        bool Saturation = false;
        bool Overlapping = false;
        unsigned int Charge = 0;
        double Path = (10.0 * thickness(DetId)) / fabs(cosine);
        double PrevGain = -1;
        double PrevGainTick = -1;
        int FirstStrip = 0;
        int NStrips = 0;

        if (StripCluster) {
          const auto& Ampls = StripCluster->amplitudes();
          FirstStrip = StripCluster->firstStrip();
          NStrips = Ampls.size();
          int APVId = FirstStrip / 128;

          if (gainHandle.isValid()) {
            PrevGain = gainHandle->getApvGain(APVId, gainHandle->getRange(DetId, 1), 1);
            PrevGainTick = gainHandle->getApvGain(APVId, gainHandle->getRange(DetId, 0), 1);
          }

          for (unsigned int a = 0; a < Ampls.size(); a++) {
            Charge += Ampls[a];
            if (Ampls[a] >= 254)
              Saturation = true;
            amplitude->push_back(Ampls[a]);
          }

          if (FirstStrip == 0)
            Overlapping = true;
          if (FirstStrip == 128)
            Overlapping = true;
          if (FirstStrip == 256)
            Overlapping = true;
          if (FirstStrip == 384)
            Overlapping = true;
          if (FirstStrip == 512)
            Overlapping = true;
          if (FirstStrip == 640)
            Overlapping = true;

          if (FirstStrip <= 127 && FirstStrip + Ampls.size() > 127)
            Overlapping = true;
          if (FirstStrip <= 255 && FirstStrip + Ampls.size() > 255)
            Overlapping = true;
          if (FirstStrip <= 383 && FirstStrip + Ampls.size() > 383)
            Overlapping = true;
          if (FirstStrip <= 511 && FirstStrip + Ampls.size() > 511)
            Overlapping = true;
          if (FirstStrip <= 639 && FirstStrip + Ampls.size() > 639)
            Overlapping = true;

          if (FirstStrip + Ampls.size() == 127)
            Overlapping = true;
          if (FirstStrip + Ampls.size() == 255)
            Overlapping = true;
          if (FirstStrip + Ampls.size() == 383)
            Overlapping = true;
          if (FirstStrip + Ampls.size() == 511)
            Overlapping = true;
          if (FirstStrip + Ampls.size() == 639)
            Overlapping = true;
          if (FirstStrip + Ampls.size() == 767)
            Overlapping = true;
        } else if (PixelCluster) {
          const auto& Ampls = PixelCluster->pixelADC();
          int FirstRow = PixelCluster->minPixelRow();
          int FirstCol = PixelCluster->minPixelCol();
          FirstStrip = ((FirstRow / 80) << 3 | (FirstCol / 52)) * 128;  //Hack to save the APVId
          NStrips = 0;
          Saturation = false;
          Overlapping = false;

          for (unsigned int a = 0; a < Ampls.size(); a++) {
            Charge += Ampls[a];
            if (Ampls[a] >= 254)
              Saturation = true;
          }
        }
#ifdef ExtendedCALIBTree
        double ChargeOverPath = (double)Charge / Path;
#endif

        trackindex->push_back(shallow::findTrackIndex(tracks, track));
        rawid->push_back(DetId);
        localdirx->push_back(trackDirection.x());
        localdiry->push_back(trackDirection.y());
        localdirz->push_back(trackDirection.z());
        firststrip->push_back(FirstStrip);
        nstrips->push_back(NStrips);
        saturation->push_back(Saturation);
        overlapping->push_back(Overlapping);
        farfromedge->push_back(StripCluster ? IsFarFromBorder(&trajState, DetId, &iSetup) : true);
        charge->push_back(Charge);
        path->push_back(Path);
#ifdef ExtendedCALIBTree
        chargeoverpath->push_back(ChargeOverPath);
#endif
        gainused->push_back(PrevGain);
        gainusedTick->push_back(PrevGainTick);
      }
    }
  }

  iEvent.put(std::move(trackindex), Prefix + "trackindex" + Suffix);
  iEvent.put(std::move(rawid), Prefix + "rawid" + Suffix);
  iEvent.put(std::move(localdirx), Prefix + "localdirx" + Suffix);
  iEvent.put(std::move(localdiry), Prefix + "localdiry" + Suffix);
  iEvent.put(std::move(localdirz), Prefix + "localdirz" + Suffix);
  iEvent.put(std::move(firststrip), Prefix + "firststrip" + Suffix);
  iEvent.put(std::move(nstrips), Prefix + "nstrips" + Suffix);
  iEvent.put(std::move(saturation), Prefix + "saturation" + Suffix);
  iEvent.put(std::move(overlapping), Prefix + "overlapping" + Suffix);
  iEvent.put(std::move(farfromedge), Prefix + "farfromedge" + Suffix);
  iEvent.put(std::move(charge), Prefix + "charge" + Suffix);
  iEvent.put(std::move(path), Prefix + "path" + Suffix);
#ifdef ExtendedCALIBTree
  iEvent.put(std::move(chargeoverpath), Prefix + "chargeoverpath" + Suffix);
#endif
  iEvent.put(std::move(amplitude), Prefix + "amplitude" + Suffix);
  iEvent.put(std::move(gainused), Prefix + "gainused" + Suffix);
  iEvent.put(std::move(gainusedTick), Prefix + "gainusedTick" + Suffix);
}

/*
void ShallowGainCalibration::beginJob(const edm::EventSetup& iSetup)
{
   printf("Befin JOB\n");

   edm::ESHandle<TrackerGeometry> tkGeom;
   iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
   vector<GeomDet*> Det = tkGeom->dets();

   edm::ESHandle<SiStripGain> gainHandle;
   iSetup.get<SiStripGainRcd>().get(gainHandle);
   if(!gainHandle.isValid()){printf("\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");exit(0);}

   for(unsigned int i=0;i<Det.size();i++){
      DetId  Detid  = Det[i]->geographicalId();
      int    SubDet = Detid.subdetId();

      if( SubDet == StripSubdetector::TIB ||  SubDet == StripSubdetector::TID ||
          SubDet == StripSubdetector::TOB ||  SubDet == StripSubdetector::TEC  ){

          StripGeomDetUnit* DetUnit     = dynamic_cast<StripGeomDetUnit*> (Det[i]);
          if(!DetUnit)continue;

          const StripTopology& Topo     = DetUnit->specificTopology();
          unsigned int         NAPV     = Topo.nstrips()/128;
   
          for(unsigned int j=0;j<NAPV;j++){
                stAPVGain* APV = new stAPVGain;
                APV->DetId         = Detid.rawId();
                APV->APVId         = j;
                APV->PreviousGain  = 1;

                APVsCollOrdered.push_back(APV);
                APVsColl[(APV->DetId<<3) | APV->APVId] = APV;
          }
      }
   }
}


void ShallowGainCalibration::beginRun(edm::Run &, const edm::EventSetup &iSetup){
    printf("BEFIN RUN\n");

    edm::ESHandle<SiStripGain> gainHandle;
    iSetup.get<SiStripGainRcd>().get(gainHandle);
    if(!gainHandle.isValid()){printf("\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");exit(0);}

    for(std::vector<stAPVGain*>::iterator it = APVsCollOrdered.begin();it!=APVsCollOrdered.end();it++){
       stAPVGain* APV = *it;
       SiStripApvGain::Range detGainRange = gainHandle->getRange(APV->DetId);
       APV->PreviousGain = *(detGainRange.first + APV->APVId);
    }
}
*/

bool ShallowGainCalibration::IsFarFromBorder(TrajectoryStateOnSurface* trajState,
                                             const uint32_t detid,
                                             const edm::EventSetup* iSetup) {
  edm::ESHandle<TrackerGeometry> tkGeom;
  iSetup->get<TrackerDigiGeometryRecord>().get(tkGeom);

  LocalPoint HitLocalPos = trajState->localPosition();
  LocalError HitLocalError = trajState->localError().positionError();

  const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
  if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
    std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
    return false;
  }

  const BoundPlane plane = it->surface();
  const TrapezoidalPlaneBounds* trapezoidalBounds(dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
  const RectangularPlaneBounds* rectangularBounds(dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));

  double DistFromBorder = 1.0;
  double HalfLength = it->surface().bounds().length() / 2.0;

  if (trapezoidalBounds) {
    std::array<const float, 4> const& parameters = (*trapezoidalBounds).parameters();
    HalfLength = parameters[3];
  } else if (rectangularBounds) {
    HalfLength = it->surface().bounds().length() / 2.0;
  } else {
    return false;
  }

  if (fabs(HitLocalPos.y()) + HitLocalError.yy() >= (HalfLength - DistFromBorder))
    return false;

  return true;
}

double ShallowGainCalibration::thickness(DetId id) {
  map<DetId, double>::iterator th = m_thicknessMap.find(id);
  if (th != m_thicknessMap.end())
    return (*th).second;
  else {
    double detThickness = 1.;
    //compute thickness normalization
    const GeomDetUnit* it = m_tracker->idToDetUnit(DetId(id));
    bool isPixel = dynamic_cast<const PixelGeomDetUnit*>(it) != nullptr;
    bool isStrip = dynamic_cast<const StripGeomDetUnit*>(it) != nullptr;
    if (!isPixel && !isStrip) {
      //FIXME throw exception
      edm::LogWarning("DeDxHitsProducer") << "\t\t this detID doesn't seem to belong to the Tracker";
      detThickness = 1.;
    } else {
      detThickness = it->surface().bounds().thickness();
    }

    m_thicknessMap[id] = detThickness;  //computed value
    return detThickness;
  }
}