DTEfficiencyTask.cc

Go to the documentation of this file.

/*
 *  See header file for a description of this class.
 *
 *  \author G. Mila - INFN Torino
 */

#include "DTEfficiencyTask.h"

// Framework
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "DQMServices/Core/interface/DQMStore.h"

//Geometry
#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"

//RecHit
#include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"

#include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
#include "CondFormats/DTObjects/interface/DTStatusFlag.h"

#include <iterator>

using namespace edm;
using namespace std;

DTEfficiencyTask::DTEfficiencyTask(const ParameterSet& pset)
    : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()), dtGeomToken_(esConsumes()) {
  debug = pset.getUntrackedParameter<bool>("debug", false);
  // the name of the 4D rec hits collection
  recHits4DToken_ =
      consumes<DTRecSegment4DCollection>(edm::InputTag(pset.getUntrackedParameter<string>("recHits4DLabel")));
  // the name of the rechits collection
  recHitToken_ = consumes<DTRecHitCollection>(edm::InputTag(pset.getUntrackedParameter<string>("recHitLabel")));

  parameters = pset;
}

DTEfficiencyTask::~DTEfficiencyTask() {}

void DTEfficiencyTask::dqmBeginRun(const edm::Run& run, const edm::EventSetup& context) {
  // Get the geometry
  muonGeom = &context.getData(muonGeomToken_);
}

void DTEfficiencyTask::bookHistograms(DQMStore::IBooker& ibooker,
                                      edm::Run const& iRun,
                                      edm::EventSetup const& context) {
  ibooker.setCurrentFolder("DT/DTEfficiencyTask");

  cout << "[DTTestPulseTask]: booking" << endl;

  //here put the static booking loop

  // Loop over all the chambers
  vector<const DTChamber*>::const_iterator ch_it = muonGeom->chambers().begin();
  vector<const DTChamber*>::const_iterator ch_end = muonGeom->chambers().end();

  for (; ch_it != ch_end; ++ch_it) {
    // Loop over the SLs
    vector<const DTSuperLayer*>::const_iterator sl_it = (*ch_it)->superLayers().begin();
    vector<const DTSuperLayer*>::const_iterator sl_end = (*ch_it)->superLayers().end();

    for (; sl_it != sl_end; ++sl_it) {
      DTSuperLayerId sl = (*sl_it)->id();
      stringstream superLayer;
      superLayer << sl.superlayer();

      // Loop over the Ls
      vector<const DTLayer*>::const_iterator l_it = (*sl_it)->layers().begin();
      vector<const DTLayer*>::const_iterator l_end = (*sl_it)->layers().end();

      for (; l_it != l_end; ++l_it) {
        DTLayerId layerId = (*l_it)->id();
        if (debug)
          cout << "   Booking histos for L: " << layerId << endl;

        // Compose the chamber name
        stringstream layer;
        layer << layerId.layer();
        stringstream superLayer;
        superLayer << layerId.superlayer();
        stringstream station;
        station << layerId.superlayerId().chamberId().station();
        stringstream sector;
        sector << layerId.superlayerId().chamberId().sector();
        stringstream wheel;
        wheel << layerId.superlayerId().chamberId().wheel();

        const int firstWire = (*l_it)->specificTopology().firstChannel();
        const int lastWire = (*l_it)->specificTopology().lastChannel();

        string lHistoName = "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str() + "_SL" +
                            superLayer.str() + "_L" + layer.str();

        ibooker.setCurrentFolder("DT/DTEfficiencyTask/Wheel" + wheel.str() + "/Station" + station.str() + "/Sector" +
                                 sector.str() + "/SuperLayer" + superLayer.str());

        // Create the monitor elements
        vector<MonitorElement*> histos;
        // histo for hits associated to the 4D reconstructed segment
        histos.push_back(ibooker.book1D("hEffOccupancy" + lHistoName,
                                        "4D segments recHits occupancy",
                                        lastWire - firstWire + 1,
                                        firstWire - 0.5,
                                        lastWire + 0.5));
        // histo for hits not associated to the segment
        histos.push_back(ibooker.book1D("hEffUnassOccupancy" + lHistoName,
                                        "4D segments recHits and Hits not associated occupancy",
                                        lastWire - firstWire + 1,
                                        firstWire - 0.5,
                                        lastWire + 0.5));
        // histo for cells associated to the 4D reconstructed segment
        histos.push_back(ibooker.book1D("hRecSegmOccupancy" + lHistoName,
                                        "4D segments cells occupancy",
                                        lastWire - firstWire + 1,
                                        firstWire - 0.5,
                                        lastWire + 0.5));

        histosPerL[layerId] = histos;

      }  // layer
    }    // superlayer
  }      // chambers
}

void DTEfficiencyTask::beginLuminosityBlock(LuminosityBlock const& lumiSeg, EventSetup const& context) {
  if (lumiSeg.id().luminosityBlock() % parameters.getUntrackedParameter<int>("ResetCycle", 3) == 0) {
    for (map<DTLayerId, vector<MonitorElement*> >::const_iterator histo = histosPerL.begin(); histo != histosPerL.end();
         histo++) {
      int size = (*histo).second.size();
      for (int i = 0; i < size; i++) {
        (*histo).second[i]->Reset();
      }
    }
  }
}

void DTEfficiencyTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
  if (debug)
    cout << "[DTEfficiencyTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event() << endl;

  // Get the 4D segment collection from the event
  edm::Handle<DTRecSegment4DCollection> all4DSegments;
  event.getByToken(recHits4DToken_, all4DSegments);

  // Get the rechit collection from the event
  Handle<DTRecHitCollection> dtRecHits;
  event.getByToken(recHitToken_, dtRecHits);

  // Get the DT Geometry
  dtGeom = &setup.getData(dtGeomToken_);

  // Loop over all chambers containing a segment
  DTRecSegment4DCollection::id_iterator chamberId;
  for (chamberId = all4DSegments->id_begin(); chamberId != all4DSegments->id_end(); ++chamberId) {
    // Get the chamber
    const DTChamber* chamber = dtGeom->chamber(*chamberId);

    // Get all 1D RecHits to be used for searches of hits not associated to segments and map them by wire
    const vector<const DTSuperLayer*>& SLayers = chamber->superLayers();
    map<DTWireId, int> wireAnd1DRecHits;
    for (vector<const DTSuperLayer*>::const_iterator superlayer = SLayers.begin(); superlayer != SLayers.end();
         superlayer++) {
      DTRecHitCollection::range range = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer((*superlayer)->id()));
      // Loop over the rechits of this ChamberId
      for (DTRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit) {
        wireAnd1DRecHits[(*rechit).wireId()] = (*rechit).wireId().wire();
      }
    }

    // Get the 4D segment range for the corresponding ChamerId
    DTRecSegment4DCollection::range range = all4DSegments->get(*chamberId);
    int nsegm = distance(range.first, range.second);
    if (debug)
      cout << "   Chamber: " << *chamberId << " has " << nsegm << " 4D segments" << endl;

    // Loop over the rechits of this ChamerId
    for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
      if (debug)
        cout << "   == RecSegment dimension: " << (*segment4D).dimension() << endl;

      // If Station != 4 skip RecHits with dimension != 4
      // For the Station 4 consider 2D RecHits
      if ((*chamberId).station() != 4 && (*segment4D).dimension() != 4) {
        if (debug)
          cout << "[DTEfficiencyTask]***Warning: RecSegment dimension is not 4 but " << (*segment4D).dimension()
               << ", skipping!" << endl;
        continue;
      } else if ((*chamberId).station() == 4 && (*segment4D).dimension() != 2) {
        if (debug)
          cout << "[DTEfficiencyTask]***Warning: RecSegment dimension is not 2 but " << (*segment4D).dimension()
               << ", skipping!" << endl;
        continue;
      }

      vector<DTRecHit1D> recHits1D;
      bool rPhi = false;
      bool rZ = false;

      // Get 1D RecHits and select only events with 7 or 8 hits in phi and 3 or 4 hits in theta (if any)
      const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
      vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
      rPhi = true;
      if (phiRecHits.size() < 7 || phiRecHits.size() > 8) {
        if (debug)
          cout << "[DTEfficiencyTask] Phi segments has: " << phiRecHits.size() << " hits, skipping"
               << endl;  // FIXME: info output
        continue;
      }
      copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D));
      const DTSLRecSegment2D* zSeg = nullptr;
      if ((*segment4D).dimension() == 4) {
        rZ = true;
        zSeg = (*segment4D).zSegment();
        vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
        if (zRecHits.size() < 3 || zRecHits.size() > 4) {
          if (debug)
            cout << "[DTEfficiencyTask] Theta segments has: " << zRecHits.size() << " hits, skipping"
                 << endl;  // FIXME: info output
          continue;
        }
        copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D));
      }

      // Skip the segment if it has more than 1 hit on the same layer
      vector<DTWireId> wireMap;
      for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D.begin(); recHit1D != recHits1D.end(); recHit1D++) {
        wireMap.push_back((*recHit1D).wireId());
      }

      bool hitsOnSameLayer = false;
      for (vector<DTWireId>::const_iterator channelId = wireMap.begin(); channelId != wireMap.end(); channelId++) {
        vector<DTWireId>::const_iterator next = channelId;
        next++;
        for (vector<DTWireId>::const_iterator ite = next; ite != wireMap.end(); ite++) {
          if ((*channelId).layerId() == (*ite).layerId()) {
            hitsOnSameLayer = true;
          }
        }
      }
      if (hitsOnSameLayer) {
        if (debug)
          cout << "[DTEfficiencyTask] This RecHit has 2 hits on the same layer, skipping!" << endl;
        continue;
      }

      // Select 2D segments with angle smaller than 45 deg
      LocalVector phiDirectionInChamber = (*phiSeg).localDirection();
      if (rPhi && fabs(phiDirectionInChamber.x() / phiDirectionInChamber.z()) > 1) {
        if (debug) {
          cout << "         RPhi segment has angle > 45 deg, skipping! " << endl;
          cout << "              Theta = " << phiDirectionInChamber.theta() << endl;
        }
        continue;
      }
      if (rZ) {
        LocalVector zDirectionInChamber = (*zSeg).localDirection();
        if (fabs(zDirectionInChamber.y() / zDirectionInChamber.z()) > 1) {
          if (debug) {
            cout << "         RZ segment has angle > 45 deg, skipping! " << endl;
            cout << "              Theta = " << zDirectionInChamber.theta() << endl;
          }
          continue;
        }
      }

      // Skip if the 4D segment has only 10 hits
      if (recHits1D.size() == 10) {
        if (debug)
          cout << "[DTEfficiencyTask] 4D Segment with only 10 hits, skipping!" << endl;
        continue;
      }

      // Analyse the case of 11 recHits for MB1,MB2,MB3 and of 7 recHits for MB4
      if ((rPhi && recHits1D.size() == 7) || (rZ && recHits1D.size() == 11)) {
        if (debug) {
          if (rPhi && recHits1D.size() == 7)
            cout << "[DTEfficiencyTask] MB4 Segment with only 7 hits!" << endl;
          if (rZ && recHits1D.size() == 11)
            cout << "[DTEfficiencyTask] 4D Segment with only 11 hits!" << endl;
        }

        // Find the layer without RecHits ----------------------------------------
        const vector<const DTSuperLayer*>& SupLayers = chamber->superLayers();
        map<DTLayerId, bool> layerMap;
        map<DTWireId, float> wireAndPosInChamberAtLayerZ;
        // Loop over layers and wires to fill a map
        for (vector<const DTSuperLayer*>::const_iterator superlayer = SupLayers.begin(); superlayer != SupLayers.end();
             superlayer++) {
          const vector<const DTLayer*> Layers = (*superlayer)->layers();
          for (vector<const DTLayer*>::const_iterator layer = Layers.begin(); layer != Layers.end(); layer++) {
            layerMap.insert(make_pair((*layer)->id(), false));
            const int firstWire = dtGeom->layer((*layer)->id())->specificTopology().firstChannel();
            const int lastWire = dtGeom->layer((*layer)->id())->specificTopology().lastChannel();
            for (int i = firstWire; i - lastWire <= 0; i++) {
              DTWireId wireId((*layer)->id(), i);
              float wireX = (*layer)->specificTopology().wirePosition(wireId.wire());
              LocalPoint wirePosInLay(wireX, 0, 0);
              GlobalPoint wirePosGlob = (*layer)->toGlobal(wirePosInLay);
              LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
              if ((*superlayer)->id().superlayer() == 1 || (*superlayer)->id().superlayer() == 3) {
                wireAndPosInChamberAtLayerZ.insert(make_pair(wireId, wirePosInChamber.x()));
              } else {
                wireAndPosInChamberAtLayerZ.insert(make_pair(wireId, wirePosInChamber.y()));
              }
            }
          }
        }

        // Loop over segment 1D RecHit
        map<DTLayerId, int> NumWireMap;
        for (vector<DTRecHit1D>::const_iterator recHit = recHits1D.begin(); recHit != recHits1D.end(); recHit++) {
          layerMap[(*recHit).wireId().layerId()] = true;
          NumWireMap[(*recHit).wireId().layerId()] = (*recHit).wireId().wire();
        }

        DTLayerId missLayerId;
        //Loop over the map and find the layer without hits
        for (map<DTLayerId, bool>::const_iterator iter = layerMap.begin(); iter != layerMap.end(); iter++) {
          if (!(*iter).second)
            missLayerId = (*iter).first;
        }
        if (debug)
          cout << "[DTEfficiencyTask] Layer without recHits is: " << missLayerId << endl;
        // -------------------------------------------------------

        const DTLayer* missLayer = chamber->layer(missLayerId);

        LocalPoint missLayerPosInChamber = chamber->toLocal(missLayer->toGlobal(LocalPoint(0, 0, 0)));

        // Segment position at Wire z in chamber local frame
        LocalPoint segPosAtZLayer = (*segment4D).localPosition() + (*segment4D).localDirection() *
                                                                       missLayerPosInChamber.z() /
                                                                       cos((*segment4D).localDirection().theta());

        DTWireId missWireId;

        // Find the id of the cell without hit ---------------------------------------------------
        for (map<DTWireId, float>::const_iterator wireAndPos = wireAndPosInChamberAtLayerZ.begin();
             wireAndPos != wireAndPosInChamberAtLayerZ.end();
             wireAndPos++) {
          DTWireId wireId = (*wireAndPos).first;
          if (wireId.layerId() == missLayerId) {
            if (missLayerId.superlayerId().superlayer() == 1 || missLayerId.superlayerId().superlayer() == 3) {
              if (fabs(segPosAtZLayer.x() - (*wireAndPos).second) < 2.1)
                missWireId = wireId;
            } else {
              if (fabs(segPosAtZLayer.y() - (*wireAndPos).second) < 2.1)
                missWireId = wireId;
            }
          }
        }
        if (debug)
          cout << "[DTEfficiencyTask] Cell without hit is: " << missWireId << endl;
        // ----------------------------------------------------------

        bool foundUnAssRechit = false;

        // Look for unassociated hits in this cell
        if (wireAnd1DRecHits.find(missWireId) != wireAnd1DRecHits.end()) {
          if (debug)
            cout << "[DTEfficiencyTask] Unassociated Hit found!" << endl;
          foundUnAssRechit = true;
        }

        for (map<DTLayerId, bool>::const_iterator iter = layerMap.begin(); iter != layerMap.end(); iter++) {
          if ((*iter).second)
            fillHistos((*iter).first,
                       dtGeom->layer((*iter).first)->specificTopology().firstChannel(),
                       dtGeom->layer((*iter).first)->specificTopology().lastChannel(),
                       NumWireMap[(*iter).first]);
          else
            fillHistos((*iter).first,
                       dtGeom->layer((*iter).first)->specificTopology().firstChannel(),
                       dtGeom->layer((*iter).first)->specificTopology().lastChannel(),
                       missWireId.wire(),
                       foundUnAssRechit);
        }

      }  // End of the loop for segment with 7 or 11 recHits

      if ((rPhi && recHits1D.size() == 8) || (rZ && recHits1D.size() == 12)) {
        map<DTLayerId, int> NumWireMap;
        DTLayerId LayerID;
        for (vector<DTRecHit1D>::const_iterator recHit = recHits1D.begin(); recHit != recHits1D.end(); recHit++) {
          LayerID = (*recHit).wireId().layerId();
          NumWireMap[LayerID] = (*recHit).wireId().wire();
        }
        for (map<DTLayerId, int>::const_iterator iter = NumWireMap.begin(); iter != NumWireMap.end(); iter++) {
          fillHistos((*iter).first,
                     dtGeom->layer((*iter).first)->specificTopology().firstChannel(),
                     dtGeom->layer((*iter).first)->specificTopology().lastChannel(),
                     NumWireMap[(*iter).first]);
        }
      }

    }  // End of loop over the 4D segments inside a sigle chamber
  }    // End of loop over all tha chamber with at least a 4D segment in the event
}

// Fill a set of histograms for a given Layer
void DTEfficiencyTask::fillHistos(DTLayerId lId, int firstWire, int lastWire, int numWire) {
  vector<MonitorElement*> histos = histosPerL[lId];
  histos[0]->Fill(numWire);
  histos[1]->Fill(numWire);
  histos[2]->Fill(numWire);
}

// Fill a set of histograms for a given Layer
void DTEfficiencyTask::fillHistos(DTLayerId lId, int firstWire, int lastWire, int missingWire, bool unassHit) {
  vector<MonitorElement*> histos = histosPerL[lId];
  if (unassHit)
    histos[1]->Fill(missingWire);
  histos[2]->Fill(missingWire);
}

// Local Variables:
// show-trailing-whitespace: t
// truncate-lines: t
// End: