DTMonitorModule/src/DTEfficiencyTask.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Mila - INFN Torino
0005  */
0006
0007 #include "DTEfficiencyTask.h"
0008
0009 // Framework
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015
0016 #include "DQMServices/Core/interface/DQMStore.h"
0017
0018 //Geometry
0019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0020 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0021
0022 //RecHit
0023 #include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
0024
0025 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
0026 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0027
0028 #include <iterator>
0029
0030 using namespace edm;
0031 using namespace std;
0032
0033 DTEfficiencyTask::DTEfficiencyTask(const ParameterSet& pset)
0034     : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()), dtGeomToken_(esConsumes()) {
0035   debug = pset.getUntrackedParameter<bool>("debug", false);
0036   // the name of the 4D rec hits collection
0037   recHits4DToken_ =
0038       consumes<DTRecSegment4DCollection>(edm::InputTag(pset.getUntrackedParameter<string>("recHits4DLabel")));
0039   // the name of the rechits collection
0040   recHitToken_ = consumes<DTRecHitCollection>(edm::InputTag(pset.getUntrackedParameter<string>("recHitLabel")));
0041
0042   parameters = pset;
0043 }
0044
0045 DTEfficiencyTask::~DTEfficiencyTask() {}
0046
0047 void DTEfficiencyTask::dqmBeginRun(const edm::Run& run, const edm::EventSetup& context) {
0048   // Get the geometry
0049   muonGeom = &context.getData(muonGeomToken_);
0050 }
0051
0052 void DTEfficiencyTask::bookHistograms(DQMStore::IBooker& ibooker,
0053                                       edm::Run const& iRun,
0054                                       edm::EventSetup const& context) {
0055   ibooker.setCurrentFolder("DT/DTEfficiencyTask");
0056
0057   cout << "[DTTestPulseTask]: booking" << endl;
0058
0059   //here put the static booking loop
0060
0061   // Loop over all the chambers
0062   vector<const DTChamber*>::const_iterator ch_it = muonGeom->chambers().begin();
0063   vector<const DTChamber*>::const_iterator ch_end = muonGeom->chambers().end();
0064
0065   for (; ch_it != ch_end; ++ch_it) {
0066     // Loop over the SLs
0067     vector<const DTSuperLayer*>::const_iterator sl_it = (*ch_it)->superLayers().begin();
0068     vector<const DTSuperLayer*>::const_iterator sl_end = (*ch_it)->superLayers().end();
0069
0070     for (; sl_it != sl_end; ++sl_it) {
0071       DTSuperLayerId sl = (*sl_it)->id();
0072       stringstream superLayer;
0073       superLayer << sl.superlayer();
0074
0075       // Loop over the Ls
0076       vector<const DTLayer*>::const_iterator l_it = (*sl_it)->layers().begin();
0077       vector<const DTLayer*>::const_iterator l_end = (*sl_it)->layers().end();
0078
0079       for (; l_it != l_end; ++l_it) {
0080         DTLayerId layerId = (*l_it)->id();
0081         if (debug)
0082           cout << "   Booking histos for L: " << layerId << endl;
0083
0084         // Compose the chamber name
0085         stringstream layer;
0086         layer << layerId.layer();
0087         stringstream superLayer;
0088         superLayer << layerId.superlayer();
0089         stringstream station;
0090         station << layerId.superlayerId().chamberId().station();
0091         stringstream sector;
0092         sector << layerId.superlayerId().chamberId().sector();
0093         stringstream wheel;
0094         wheel << layerId.superlayerId().chamberId().wheel();
0095
0096         const int firstWire = (*l_it)->specificTopology().firstChannel();
0097         const int lastWire = (*l_it)->specificTopology().lastChannel();
0098
0099         string lHistoName = "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str() + "_SL" +
0100                             superLayer.str() + "_L" + layer.str();
0101
0102         ibooker.setCurrentFolder("DT/DTEfficiencyTask/Wheel" + wheel.str() + "/Station" + station.str() + "/Sector" +
0103                                  sector.str() + "/SuperLayer" + superLayer.str());
0104
0105         // Create the monitor elements
0106         vector<MonitorElement*> histos;
0107         // histo for hits associated to the 4D reconstructed segment
0108         histos.push_back(ibooker.book1D("hEffOccupancy" + lHistoName,
0109                                         "4D segments recHits occupancy",
0110                                         lastWire - firstWire + 1,
0111                                         firstWire - 0.5,
0112                                         lastWire + 0.5));
0113         // histo for hits not associated to the segment
0114         histos.push_back(ibooker.book1D("hEffUnassOccupancy" + lHistoName,
0115                                         "4D segments recHits and Hits not associated occupancy",
0116                                         lastWire - firstWire + 1,
0117                                         firstWire - 0.5,
0118                                         lastWire + 0.5));
0119         // histo for cells associated to the 4D reconstructed segment
0120         histos.push_back(ibooker.book1D("hRecSegmOccupancy" + lHistoName,
0121                                         "4D segments cells occupancy",
0122                                         lastWire - firstWire + 1,
0123                                         firstWire - 0.5,
0124                                         lastWire + 0.5));
0125
0126         histosPerL[layerId] = histos;
0127
0128       }  // layer
0129     }  // superlayer
0130   }  // chambers
0131 }
0132
0133 void DTEfficiencyTask::beginLuminosityBlock(LuminosityBlock const& lumiSeg, EventSetup const& context) {
0134   if (lumiSeg.id().luminosityBlock() % parameters.getUntrackedParameter<int>("ResetCycle", 3) == 0) {
0135     for (map<DTLayerId, vector<MonitorElement*> >::const_iterator histo = histosPerL.begin(); histo != histosPerL.end();
0136          histo++) {
0137       int size = (*histo).second.size();
0138       for (int i = 0; i < size; i++) {
0139         (*histo).second[i]->Reset();
0140       }
0141     }
0142   }
0143 }
0144
0145 void DTEfficiencyTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0146   if (debug)
0147     cout << "[DTEfficiencyTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event() << endl;
0148
0149   // Get the 4D segment collection from the event
0150   edm::Handle<DTRecSegment4DCollection> all4DSegments;
0151   event.getByToken(recHits4DToken_, all4DSegments);
0152
0153   // Get the rechit collection from the event
0154   Handle<DTRecHitCollection> dtRecHits;
0155   event.getByToken(recHitToken_, dtRecHits);
0156
0157   // Get the DT Geometry
0158   dtGeom = &setup.getData(dtGeomToken_);
0159
0160   // Loop over all chambers containing a segment
0161   DTRecSegment4DCollection::id_iterator chamberId;
0162   for (chamberId = all4DSegments->id_begin(); chamberId != all4DSegments->id_end(); ++chamberId) {
0163     // Get the chamber
0164     const DTChamber* chamber = dtGeom->chamber(*chamberId);
0165
0166     // Get all 1D RecHits to be used for searches of hits not associated to segments and map them by wire
0167     const vector<const DTSuperLayer*>& SLayers = chamber->superLayers();
0168     map<DTWireId, int> wireAnd1DRecHits;
0169     for (vector<const DTSuperLayer*>::const_iterator superlayer = SLayers.begin(); superlayer != SLayers.end();
0170          superlayer++) {
0171       DTRecHitCollection::range range = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer((*superlayer)->id()));
0172       // Loop over the rechits of this ChamberId
0173       for (DTRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit) {
0174         wireAnd1DRecHits[(*rechit).wireId()] = (*rechit).wireId().wire();
0175       }
0176     }
0177
0178     // Get the 4D segment range for the corresponding ChamerId
0179     DTRecSegment4DCollection::range range = all4DSegments->get(*chamberId);
0180     int nsegm = distance(range.first, range.second);
0181     if (debug)
0182       cout << "   Chamber: " << *chamberId << " has " << nsegm << " 4D segments" << endl;
0183
0184     // Loop over the rechits of this ChamerId
0185     for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0186       if (debug)
0187         cout << "   == RecSegment dimension: " << (*segment4D).dimension() << endl;
0188
0189       // If Station != 4 skip RecHits with dimension != 4
0190       // For the Station 4 consider 2D RecHits
0191       if ((*chamberId).station() != 4 && (*segment4D).dimension() != 4) {
0192         if (debug)
0193           cout << "[DTEfficiencyTask]***Warning: RecSegment dimension is not 4 but " << (*segment4D).dimension()
0194                << ", skipping!" << endl;
0195         continue;
0196       } else if ((*chamberId).station() == 4 && (*segment4D).dimension() != 2) {
0197         if (debug)
0198           cout << "[DTEfficiencyTask]***Warning: RecSegment dimension is not 2 but " << (*segment4D).dimension()
0199                << ", skipping!" << endl;
0200         continue;
0201       }
0202
0203       vector<DTRecHit1D> recHits1D;
0204       bool rPhi = false;
0205       bool rZ = false;
0206
0207       // Get 1D RecHits and select only events with 7 or 8 hits in phi and 3 or 4 hits in theta (if any)
0208       const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
0209       vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
0210       rPhi = true;
0211       if (phiRecHits.size() < 7 || phiRecHits.size() > 8) {
0212         if (debug)
0213           cout << "[DTEfficiencyTask] Phi segments has: " << phiRecHits.size() << " hits, skipping"
0214                << endl;  // FIXME: info output
0215         continue;
0216       }
0217       copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D));
0218       const DTSLRecSegment2D* zSeg = nullptr;
0219       if ((*segment4D).dimension() == 4) {
0220         rZ = true;
0221         zSeg = (*segment4D).zSegment();
0222         vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
0223         if (zRecHits.size() < 3 || zRecHits.size() > 4) {
0224           if (debug)
0225             cout << "[DTEfficiencyTask] Theta segments has: " << zRecHits.size() << " hits, skipping"
0226                  << endl;  // FIXME: info output
0227           continue;
0228         }
0229         copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D));
0230       }
0231
0232       // Skip the segment if it has more than 1 hit on the same layer
0233       vector<DTWireId> wireMap;
0234       for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D.begin(); recHit1D != recHits1D.end(); recHit1D++) {
0235         wireMap.push_back((*recHit1D).wireId());
0236       }
0237
0238       bool hitsOnSameLayer = false;
0239       for (vector<DTWireId>::const_iterator channelId = wireMap.begin(); channelId != wireMap.end(); channelId++) {
0240         vector<DTWireId>::const_iterator next = channelId;
0241         next++;
0242         for (vector<DTWireId>::const_iterator ite = next; ite != wireMap.end(); ite++) {
0243           if ((*channelId).layerId() == (*ite).layerId()) {
0244             hitsOnSameLayer = true;
0245           }
0246         }
0247       }
0248       if (hitsOnSameLayer) {
0249         if (debug)
0250           cout << "[DTEfficiencyTask] This RecHit has 2 hits on the same layer, skipping!" << endl;
0251         continue;
0252       }
0253
0254       // Select 2D segments with angle smaller than 45 deg
0255       LocalVector phiDirectionInChamber = (*phiSeg).localDirection();
0256       if (rPhi && fabs(phiDirectionInChamber.x() / phiDirectionInChamber.z()) > 1) {
0257         if (debug) {
0258           cout << "         RPhi segment has angle > 45 deg, skipping! " << endl;
0259           cout << "              Theta = " << phiDirectionInChamber.theta() << endl;
0260         }
0261         continue;
0262       }
0263       if (rZ) {
0264         LocalVector zDirectionInChamber = (*zSeg).localDirection();
0265         if (fabs(zDirectionInChamber.y() / zDirectionInChamber.z()) > 1) {
0266           if (debug) {
0267             cout << "         RZ segment has angle > 45 deg, skipping! " << endl;
0268             cout << "              Theta = " << zDirectionInChamber.theta() << endl;
0269           }
0270           continue;
0271         }
0272       }
0273
0274       // Skip if the 4D segment has only 10 hits
0275       if (recHits1D.size() == 10) {
0276         if (debug)
0277           cout << "[DTEfficiencyTask] 4D Segment with only 10 hits, skipping!" << endl;
0278         continue;
0279       }
0280
0281       // Analyse the case of 11 recHits for MB1,MB2,MB3 and of 7 recHits for MB4
0282       if ((rPhi && recHits1D.size() == 7) || (rZ && recHits1D.size() == 11)) {
0283         if (debug) {
0284           if (rPhi && recHits1D.size() == 7)
0285             cout << "[DTEfficiencyTask] MB4 Segment with only 7 hits!" << endl;
0286           if (rZ && recHits1D.size() == 11)
0287             cout << "[DTEfficiencyTask] 4D Segment with only 11 hits!" << endl;
0288         }
0289
0290         // Find the layer without RecHits ----------------------------------------
0291         const vector<const DTSuperLayer*>& SupLayers = chamber->superLayers();
0292         map<DTLayerId, bool> layerMap;
0293         map<DTWireId, float> wireAndPosInChamberAtLayerZ;
0294         // Loop over layers and wires to fill a map
0295         for (vector<const DTSuperLayer*>::const_iterator superlayer = SupLayers.begin(); superlayer != SupLayers.end();
0296              superlayer++) {
0297           const vector<const DTLayer*> Layers = (*superlayer)->layers();
0298           for (vector<const DTLayer*>::const_iterator layer = Layers.begin(); layer != Layers.end(); layer++) {
0299             layerMap.insert(make_pair((*layer)->id(), false));
0300             const int firstWire = dtGeom->layer((*layer)->id())->specificTopology().firstChannel();
0301             const int lastWire = dtGeom->layer((*layer)->id())->specificTopology().lastChannel();
0302             for (int i = firstWire; i - lastWire <= 0; i++) {
0303               DTWireId wireId((*layer)->id(), i);
0304               float wireX = (*layer)->specificTopology().wirePosition(wireId.wire());
0305               LocalPoint wirePosInLay(wireX, 0, 0);
0306               GlobalPoint wirePosGlob = (*layer)->toGlobal(wirePosInLay);
0307               LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
0308               if ((*superlayer)->id().superlayer() == 1 || (*superlayer)->id().superlayer() == 3) {
0309                 wireAndPosInChamberAtLayerZ.insert(make_pair(wireId, wirePosInChamber.x()));
0310               } else {
0311                 wireAndPosInChamberAtLayerZ.insert(make_pair(wireId, wirePosInChamber.y()));
0312               }
0313             }
0314           }
0315         }
0316
0317         // Loop over segment 1D RecHit
0318         map<DTLayerId, int> NumWireMap;
0319         for (vector<DTRecHit1D>::const_iterator recHit = recHits1D.begin(); recHit != recHits1D.end(); recHit++) {
0320           layerMap[(*recHit).wireId().layerId()] = true;
0321           NumWireMap[(*recHit).wireId().layerId()] = (*recHit).wireId().wire();
0322         }
0323
0324         DTLayerId missLayerId;
0325         //Loop over the map and find the layer without hits
0326         for (map<DTLayerId, bool>::const_iterator iter = layerMap.begin(); iter != layerMap.end(); iter++) {
0327           if (!(*iter).second)
0328             missLayerId = (*iter).first;
0329         }
0330         if (debug)
0331           cout << "[DTEfficiencyTask] Layer without recHits is: " << missLayerId << endl;
0332         // -------------------------------------------------------
0333
0334         const DTLayer* missLayer = chamber->layer(missLayerId);
0335
0336         LocalPoint missLayerPosInChamber = chamber->toLocal(missLayer->toGlobal(LocalPoint(0, 0, 0)));
0337
0338         // Segment position at Wire z in chamber local frame
0339         LocalPoint segPosAtZLayer = (*segment4D).localPosition() + (*segment4D).localDirection() *
0340                                                                        missLayerPosInChamber.z() /
0341                                                                        cos((*segment4D).localDirection().theta());
0342
0343         DTWireId missWireId;
0344
0345         // Find the id of the cell without hit ---------------------------------------------------
0346         for (map<DTWireId, float>::const_iterator wireAndPos = wireAndPosInChamberAtLayerZ.begin();
0347              wireAndPos != wireAndPosInChamberAtLayerZ.end();
0348              wireAndPos++) {
0349           DTWireId wireId = (*wireAndPos).first;
0350           if (wireId.layerId() == missLayerId) {
0351             if (missLayerId.superlayerId().superlayer() == 1 || missLayerId.superlayerId().superlayer() == 3) {
0352               if (fabs(segPosAtZLayer.x() - (*wireAndPos).second) < 2.1)
0353                 missWireId = wireId;
0354             } else {
0355               if (fabs(segPosAtZLayer.y() - (*wireAndPos).second) < 2.1)
0356                 missWireId = wireId;
0357             }
0358           }
0359         }
0360         if (debug)
0361           cout << "[DTEfficiencyTask] Cell without hit is: " << missWireId << endl;
0362         // ----------------------------------------------------------
0363
0364         bool foundUnAssRechit = false;
0365
0366         // Look for unassociated hits in this cell
0367         if (wireAnd1DRecHits.find(missWireId) != wireAnd1DRecHits.end()) {
0368           if (debug)
0369             cout << "[DTEfficiencyTask] Unassociated Hit found!" << endl;
0370           foundUnAssRechit = true;
0371         }
0372
0373         for (map<DTLayerId, bool>::const_iterator iter = layerMap.begin(); iter != layerMap.end(); iter++) {
0374           if ((*iter).second)
0375             fillHistos((*iter).first,
0376                        dtGeom->layer((*iter).first)->specificTopology().firstChannel(),
0377                        dtGeom->layer((*iter).first)->specificTopology().lastChannel(),
0378                        NumWireMap[(*iter).first]);
0379           else
0380             fillHistos((*iter).first,
0381                        dtGeom->layer((*iter).first)->specificTopology().firstChannel(),
0382                        dtGeom->layer((*iter).first)->specificTopology().lastChannel(),
0383                        missWireId.wire(),
0384                        foundUnAssRechit);
0385         }
0386
0387       }  // End of the loop for segment with 7 or 11 recHits
0388
0389       if ((rPhi && recHits1D.size() == 8) || (rZ && recHits1D.size() == 12)) {
0390         map<DTLayerId, int> NumWireMap;
0391         DTLayerId LayerID;
0392         for (vector<DTRecHit1D>::const_iterator recHit = recHits1D.begin(); recHit != recHits1D.end(); recHit++) {
0393           LayerID = (*recHit).wireId().layerId();
0394           NumWireMap[LayerID] = (*recHit).wireId().wire();
0395         }
0396         for (map<DTLayerId, int>::const_iterator iter = NumWireMap.begin(); iter != NumWireMap.end(); iter++) {
0397           fillHistos((*iter).first,
0398                      dtGeom->layer((*iter).first)->specificTopology().firstChannel(),
0399                      dtGeom->layer((*iter).first)->specificTopology().lastChannel(),
0400                      NumWireMap[(*iter).first]);
0401         }
0402       }
0403
0404     }  // End of loop over the 4D segments inside a sigle chamber
0405   }  // End of loop over all tha chamber with at least a 4D segment in the event
0406 }
0407
0408 // Fill a set of histograms for a given Layer
0409 void DTEfficiencyTask::fillHistos(DTLayerId lId, int firstWire, int lastWire, int numWire) {
0410   vector<MonitorElement*> histos = histosPerL[lId];
0411   histos[0]->Fill(numWire);
0412   histos[1]->Fill(numWire);
0413   histos[2]->Fill(numWire);
0414 }
0415
0416 // Fill a set of histograms for a given Layer
0417 void DTEfficiencyTask::fillHistos(DTLayerId lId, int firstWire, int lastWire, int missingWire, bool unassHit) {
0418   vector<MonitorElement*> histos = histosPerL[lId];
0419   if (unassHit)
0420     histos[1]->Fill(missingWire);
0421   histos[2]->Fill(missingWire);
0422 }
0423
0424 // Local Variables:
0425 // show-trailing-whitespace: t
0426 // truncate-lines: t
0427 // End: