DTMonitorModule/src/DTChamberEfficiency.cc

0001 /******* \class DTEffAnalyzer *******
0002  *
0003  * Description:
0004  *
0005  *  detailed description
0006  *
0007  * \author : Mario Pelliccioni, pellicci@cern.ch
0008  * $date   : 20/11/2008 16:50:57 CET $
0009  *
0010  * Modification:
0011  *
0012  *********************************/
0013
0014 #include "DTChamberEfficiency.h"
0015
0016 #include "DataFormats/Common/interface/Handle.h"
0017
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021
0022 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
0023 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0024
0025 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0026 #include "MagneticField/Engine/interface/MagneticField.h"
0027
0028 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
0029
0030 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0031
0032 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0033
0034 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0035 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0036 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0037 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0038
0039 #include "DataFormats/DetId/interface/DetId.h"
0040 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0041 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0042 #include "DataFormats/Common/interface/RefToBase.h"
0043
0044 #include <cmath>
0045
0046 using namespace std;
0047 using namespace edm;
0048
0049 DTChamberEfficiency::DTChamberEfficiency(const ParameterSet& pSet) {
0050   // Get the debug parameter for verbose output
0051   debug = pSet.getUntrackedParameter<bool>("debug", false);
0052
0053   LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiency") << "DTChamberEfficiency: constructor called";
0054
0055   // service parameters
0056   ParameterSet serviceParameters = pSet.getParameter<ParameterSet>("ServiceParameters");
0057   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0058
0059   theTracksLabel_ = pSet.getUntrackedParameter<InputTag>("TrackCollection");
0060   theTracksToken_ = consumes<reco::TrackCollection>(theTracksLabel_);
0061
0062   theMaxChi2 = static_cast<unsigned int>(pSet.getParameter<double>("theMaxChi2"));
0063   theNSigma = pSet.getParameter<double>("theNSigma");
0064   theMinNrec = static_cast<int>(pSet.getParameter<double>("theMinNrec"));
0065
0066   labelRPCRecHits = pSet.getUntrackedParameter<InputTag>("theRPCRecHits");
0067
0068   thedt4DSegments = pSet.getUntrackedParameter<InputTag>("dt4DSegments");
0069   thecscSegments = pSet.getUntrackedParameter<InputTag>("cscSegments");
0070
0071   edm::ConsumesCollector iC = consumesCollector();
0072
0073   theMeasurementExtractor = new MuonDetLayerMeasurements(
0074       thedt4DSegments, thecscSegments, labelRPCRecHits, InputTag(), InputTag(), iC, true, false, false, false);
0075
0076   theNavigationType = pSet.getParameter<string>("NavigationType");
0077
0078   theEstimator = new Chi2MeasurementEstimator(theMaxChi2, theNSigma);
0079 }
0080
0081 DTChamberEfficiency::~DTChamberEfficiency() {
0082   LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") << "DTChamberEfficiency: destructor called";
0083
0084   // free memory
0085   delete theService;
0086   delete theMeasurementExtractor;
0087   delete theEstimator;
0088 }
0089
0090 void DTChamberEfficiency::dqmBeginRun(const Run&, const EventSetup&) {}
0091
0092 void DTChamberEfficiency::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const&) {
0093   LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") << "DTChamberEfficiency: booking histos";
0094
0095   // Create the monitor elements
0096   ibooker.setCurrentFolder("DT/05-ChamberEff/Task");
0097
0098   for (int wheel = -2; wheel <= 2; wheel++) {
0099     vector<MonitorElement*> histos;
0100
0101     stringstream wheel_str;
0102     wheel_str << wheel;
0103
0104     histos.push_back(ibooker.book2D(
0105         "hCountSectVsChamb_All_W" + wheel_str.str(), "Countings for wheel " + wheel_str.str(), 14, 1., 15., 4, 1., 5.));
0106
0107     histos.push_back(ibooker.book2D(
0108         "hCountSectVsChamb_Qual_W" + wheel_str.str(), "Countings for wheel " + wheel_str.str(), 14, 1., 15., 4, 1., 5.));
0109
0110     histos.push_back(ibooker.book2D(
0111         "hExtrapSectVsChamb_W" + wheel_str.str(), "Extrapolations for wheel " + wheel_str.str(), 14, 1., 15., 4, 1., 5.));
0112
0113     histosPerW.push_back(histos);
0114   }
0115
0116   return;
0117 }
0118
0119 void DTChamberEfficiency::analyze(const Event& event, const EventSetup& eventSetup) {
0120   LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency")
0121       << "--- [DTChamberEfficiency] Event analysed #Run: " << event.id().run() << " #Event: " << event.id().event()
0122       << endl;
0123
0124   theService->update(eventSetup);
0125   theMeasurementExtractor->setEvent(event);
0126
0127   //Read tracks from event
0128   Handle<reco::TrackCollection> tracks;
0129   event.getByToken(theTracksToken_, tracks);
0130
0131   if (tracks.isValid()) {  // check the validity of the collection
0132
0133     const edm::ESHandle<GlobalTrackingGeometry>& globalTrackingGeometry = theService->trackingGeometry();
0134     const MagneticField* magneticField = theService->magneticField().product();
0135
0136     //loop over the muons
0137     for (reco::TrackCollection::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
0138       reco::TransientTrack trans_track(*track, magneticField, globalTrackingGeometry);
0139       const int recHitsize = (int)trans_track.recHitsSize();
0140       if (recHitsize < theMinNrec)
0141         continue;
0142
0143       // printout the DT rechits used by the track
0144       if (debug) {
0145         LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") << "--- New track" << endl;
0146         set<DTChamberId> chAlrUsed;
0147         for (trackingRecHit_iterator rHit = trans_track.recHitsBegin(); rHit != trans_track.recHitsEnd(); ++rHit) {
0148           DetId rHitid = (*rHit)->geographicalId();
0149           if (!(rHitid.det() == DetId::Muon && rHitid.subdetId() == MuonSubdetId::DT))
0150             continue;
0151           DTChamberId wId(rHitid.rawId());
0152           if (chAlrUsed.find(wId) != chAlrUsed.end())
0153             continue;
0154           chAlrUsed.insert(wId);
0155           LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") << "   " << wId << endl;
0156         }
0157       }
0158
0159       // Get the layer on which the seed relies
0160       DetId id = trans_track.track().innerDetId();
0161       const DetLayer* initialLayer = theService->detLayerGeometry()->idToLayer(id);
0162
0163       TrajectoryStateOnSurface init_fs = trans_track.innermostMeasurementState();
0164       const FreeTrajectoryState* init_fs_free = init_fs.freeState();
0165
0166       //get the list of compatible layers
0167       vector<const DetLayer*> layer_list =
0168           compatibleLayers(*theService->muonNavigationSchool(), initialLayer, *init_fs_free, alongMomentum);
0169       vector<const DetLayer*> layer_list_2 =
0170           compatibleLayers(*theService->muonNavigationSchool(), initialLayer, *init_fs_free, oppositeToMomentum);
0171
0172       layer_list.insert(layer_list.end(), layer_list_2.begin(), layer_list_2.end());
0173
0174       set<DTChamberId> alreadyCheckedCh;
0175
0176       //loop over the list of compatible layers
0177       for (int i = 0; i < (int)layer_list.size(); i++) {
0178         //propagate the track to the i-th layer
0179         TrajectoryStateOnSurface tsos = propagator()->propagate(init_fs, layer_list.at(i)->surface());
0180         if (!tsos.isValid())
0181           continue;
0182
0183         //determine the chambers kinematically compatible with the track on the i-th layer
0184         vector<DetWithState> dss = layer_list.at(i)->compatibleDets(tsos, *propagator(), *theEstimator);
0185
0186         if (dss.empty())
0187           continue;
0188
0189         // get the first det (it's the most compatible)
0190         const DetWithState detWithState = dss.front();
0191         const DetId idDetLay = detWithState.first->geographicalId();
0192
0193         // check if this is a DT and the track has the needed quality
0194         if (!chamberSelection(idDetLay, trans_track))
0195           continue;
0196
0197         DTChamberId DTid = (DTChamberId)idDetLay;
0198
0199         // check if the chamber has already been counted
0200         if (alreadyCheckedCh.find(DTid) != alreadyCheckedCh.end())
0201           continue;
0202         alreadyCheckedCh.insert(DTid);
0203
0204         // get the compatible measurements
0205         MeasurementContainer detMeasurements_initial = theMeasurementExtractor->measurements(
0206             layer_list.at(i), detWithState.first, detWithState.second, *theEstimator, event);
0207         LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency")
0208             << "     chamber: " << DTid << " has: " << detMeasurements_initial.size() << " comp. meas." << endl;
0209
0210         //we want to be more picky about the quality of the segments:
0211         //exclude the segments with less than 12 hits
0212         MeasurementContainer detMeasurements = segQualityCut(detMeasurements_initial);
0213
0214         // get the histos for this chamber
0215         vector<MonitorElement*> histos = histosPerW[DTid.wheel() + 2];
0216         // fill them
0217         if (!detMeasurements_initial.empty())
0218           histos[0]->Fill(DTid.sector(), DTid.station(), 1.);
0219         if (!detMeasurements.empty())
0220           histos[1]->Fill(DTid.sector(), DTid.station(), 1.);
0221         histos[2]->Fill(DTid.sector(), DTid.station(), 1.);
0222       }
0223     }
0224   } else {
0225     LogInfo("DTDQM|DTMonitorModule|DTChamberEfficiency")
0226         << "[DTChamberEfficiency] Collection: " << theTracksLabel_ << " is not valid!" << endl;
0227   }
0228   return;
0229 }
0230
0231 bool DTChamberEfficiency::chamberSelection(const DetId& idDetLay, reco::TransientTrack& trans_track) const {
0232   //check that we have a muon and that is a DT detector
0233   if (!(idDetLay.det() == DetId::Muon && idDetLay.subdetId() == MuonSubdetId::DT))
0234     return false;
0235
0236   if (trans_track.recHitsSize() == 2)
0237     if (trans_track.recHit(0)->geographicalId() == idDetLay || trans_track.recHit(1)->geographicalId() == idDetLay)
0238       return false;
0239
0240   return true;
0241 }
0242
0243 MeasurementContainer DTChamberEfficiency::segQualityCut(const MeasurementContainer& seg_list) const {
0244   MeasurementContainer result;
0245
0246   for (MeasurementContainer::const_iterator mescont_Itr = seg_list.begin(); mescont_Itr != seg_list.end();
0247        ++mescont_Itr) {
0248     //get the rechits of the segment
0249     TransientTrackingRecHit::ConstRecHitContainer recHit_list = mescont_Itr->recHit()->transientHits();
0250
0251     //loop over the rechits and get the number of hits
0252     int nhit_seg(0);
0253     for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator recList_Itr = recHit_list.begin();
0254          recList_Itr != recHit_list.end();
0255          ++recList_Itr) {
0256       nhit_seg += (int)(*recList_Itr)->transientHits().size();
0257     }
0258
0259     DTChamberId tmpId = (DTChamberId)mescont_Itr->recHit()->hit()->geographicalId();
0260
0261     if (tmpId.station() < 4 && nhit_seg >= 12)
0262       result.push_back(*mescont_Itr);
0263     if (tmpId.station() == 4 && nhit_seg >= 8)
0264       result.push_back(*mescont_Itr);
0265   }
0266
0267   return result;
0268 }
0269
0270 vector<const DetLayer*> DTChamberEfficiency::compatibleLayers(const NavigationSchool& navigationSchool,
0271                                                               const DetLayer* initialLayer,
0272                                                               const FreeTrajectoryState& fts,
0273                                                               PropagationDirection propDir) {
0274   vector<const DetLayer*> detLayers;
0275
0276   if (theNavigationType == "Standard") {
0277     // ask for compatible layers
0278     detLayers = navigationSchool.compatibleLayers(*initialLayer, fts, propDir);
0279     // I have to fit by hand the first layer until the seedTSOS is defined on the first rechit layer
0280     // In fact the first layer is not returned by initialLayer->compatibleLayers.
0281
0282     detLayers.insert(detLayers.begin(), initialLayer);
0283
0284   } else if (theNavigationType == "Direct") {
0285     DirectMuonNavigation navigation(ESHandle<MuonDetLayerGeometry>(&*theService->detLayerGeometry()));
0286     detLayers = navigation.compatibleLayers(fts, propDir);
0287   } else
0288     LogError("DTDQM|DTMonitorModule|DTChamberEfficiency") << "No Properly Navigation Selected!!" << endl;
0289
0290   return detLayers;
0291 }
0292
0293 inline ESHandle<Propagator> DTChamberEfficiency::propagator() const {
0294   return theService->propagator("SteppingHelixPropagatorAny");
0295 }
0296
0297 // Local Variables:
0298 // show-trailing-whitespace: t
0299 // truncate-lines: t
0300 // End: