Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-21 22:39:42

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