CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DTChamberEfficiency.cc
Go to the documentation of this file.
1 /******* \class DTEffAnalyzer *******
2  *
3  * Description:
4  *
5  * detailed description
6  *
7  * \author : Mario Pelliccioni, pellicci@cern.ch
8  * $date : 20/11/2008 16:50:57 CET $
9  *
10  * Modification:
11  *
12  *********************************/
13 
14 #include "DTChamberEfficiency.h"
15 
20 
22 
24 
27 
30 
32 
34 
36 
38 
39 
44 
49 
50 
51 
52 
53 
54 #include <cmath>
55 
56 using namespace std;
57 using namespace edm;
58 
60 {
61  // Get the debug parameter for verbose output
62  debug = pSet.getUntrackedParameter<bool>("debug","false");
63 
64  LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiency")
65  << "DTChamberEfficiency: constructor called";
66 
67  // Get the DQM needed services
68  theDbe = Service<DQMStore>().operator->();
69  theDbe->setCurrentFolder("DT/05-ChamberEff/Task");
70 
71  // service parameters
72  ParameterSet serviceParameters = pSet.getParameter<ParameterSet>("ServiceParameters");
73  theService = new MuonServiceProxy(serviceParameters);
74 
75  theTracksLabel = pSet.getParameter<InputTag>("TrackCollection");
76 
77  theMaxChi2 = static_cast<unsigned int>(pSet.getParameter<double>("theMaxChi2"));
78  theNSigma = pSet.getParameter<double>("theNSigma");
79  theMinNrec = static_cast<int>(pSet.getParameter<double>("theMinNrec"));
80 
81  labelRPCRecHits = pSet.getParameter<InputTag>("theRPCRecHits");
82 
83  thedt4DSegments = pSet.getParameter<InputTag>("dt4DSegments");
84  thecscSegments = pSet.getParameter<InputTag>("cscSegments");
85 
86  theMeasurementExtractor = new MuonDetLayerMeasurements(thedt4DSegments,thecscSegments,
87  labelRPCRecHits,true,false,false);
88 
89  theNavigationType = pSet.getParameter<string>("NavigationType");
90 
91  theEstimator = new Chi2MeasurementEstimator(theMaxChi2,theNSigma);
92 }
93 
94 
96 {
97  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency")
98  << "DTChamberEfficiency: destructor called";
99 
100  // free memory
101  delete theService;
102  delete theMeasurementExtractor;
103  delete theEstimator;
104 }
105 
107 
108  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency")
109  << "DTChamberEfficiency: beginOfJob";
110 
111  bookHistos();
112 
113  return;
114 }
115 
117 {
118  // Get the DT Geometry
119  setup.get<MuonGeometryRecord>().get(dtGeom);
120 
121  setup.get<IdealMagneticFieldRecord>().get(magfield);
122  setup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
123 
124  return;
125 }
126 
128 {
129  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency")
130  << "DTChamberEfficiency: endOfJob";
131 
132  return;
133 }
134 
135 // Book a set of histograms for a given Layer
137 {
138  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency")
139  << "DTChamberEfficiency: booking histos";
140 
141  // Create the monitor elements
142  theDbe->setCurrentFolder("DT/05-ChamberEff/Task");
143 
144  for(int wheel=-2;wheel<=2;wheel++){
145 
146  vector<MonitorElement *> histos;
147 
148  stringstream wheel_str; wheel_str << wheel;
149 
150  histos.push_back(theDbe->book2D("hCountSectVsChamb_All_W"+ wheel_str.str(),
151  "Countings for wheel " + wheel_str.str(),14,1.,15.,4,1.,5.));
152 
153  histos.push_back(theDbe->book2D("hCountSectVsChamb_Qual_W"+ wheel_str.str(),
154  "Countings for wheel " + wheel_str.str(),14,1.,15.,4,1.,5.));
155 
156 
157  histos.push_back(theDbe->book2D("hExtrapSectVsChamb_W"+ wheel_str.str(),
158  "Extrapolations for wheel " + wheel_str.str(),14,1.,15.,4,1.,5.));
159 
160  histosPerW.push_back(histos);
161  }
162 
163  return;
164 }
165 
167  const EventSetup& eventSetup)
168 {
169 
170  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") <<
171  "--- [DTChamberEfficiency] Event analysed #Run: " <<
172  event.id().run() << " #Event: " << event.id().event() << endl;
173 
174  theService->update(eventSetup);
175  theMeasurementExtractor->setEvent(event);
176 
177  //Read tracks from event
179  event.getByLabel(theTracksLabel, tracks);
180 
181  if(tracks.isValid()) { // check the validity of the collection
182 
183  //loop over the muons
184  for(reco::TrackCollection::const_iterator track = tracks->begin(); track!=tracks->end(); ++track) {
185 
186 
187  reco::TransientTrack trans_track(*track,magfield.product(),theTrackingGeometry);
188  const int recHitsize = (int)trans_track.recHitsSize();
189  if(recHitsize < theMinNrec) continue;
190 
191  // printout the DT rechits used by the track
192  if (debug) {
193  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") << "--- New track" << endl;
194  set<DTChamberId> chAlrUsed;
195  for(trackingRecHit_iterator rHit = trans_track.recHitsBegin();
196  rHit != trans_track.recHitsEnd(); ++rHit) {
197  DetId rHitid = (*rHit)->geographicalId();
198  if(!(rHitid.det() == DetId::Muon && rHitid.subdetId() == MuonSubdetId::DT)) continue;
199  DTChamberId wId(rHitid.rawId());
200  if(chAlrUsed.find(wId) != chAlrUsed.end()) continue;
201  chAlrUsed.insert(wId);
202  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") << " " << wId << endl;
203  }
204  }
205 
206  // Get the layer on which the seed relies
207  DetId id = trans_track.track().innerDetId();
208  const DetLayer *initialLayer = theService->detLayerGeometry()->idToLayer(id);
209 
210  TrajectoryStateOnSurface init_fs = trans_track.innermostMeasurementState();
211  FreeTrajectoryState *init_fs_free = init_fs.freeState();
212 
213  //get the list of compatible layers
214  vector<const DetLayer*> layer_list = compatibleLayers(initialLayer,*init_fs_free,alongMomentum);
215  vector<const DetLayer*> layer_list_2 = compatibleLayers(initialLayer,*init_fs_free,oppositeToMomentum);
216 
217  layer_list.insert(layer_list.end(),layer_list_2.begin(),layer_list_2.end());
218 
219  set<DTChamberId> alreadyCheckedCh;
220 
221  //loop over the list of compatible layers
222  for(int i=0;i< (int)layer_list.size();i++) {
223 
224  //propagate the track to the i-th layer
225  TrajectoryStateOnSurface tsos = propagator()->propagate(init_fs,layer_list.at(i)->surface());
226  if(!tsos.isValid()) continue;
227 
228  //determine the chambers kinematically compatible with the track on the i-th layer
229  vector<DetWithState> dss = layer_list.at(i)->compatibleDets(tsos, *propagator(), *theEstimator);
230 
231  if(dss.size() == 0) continue;
232 
233  // get the first det (it's the most compatible)
234  const DetWithState detWithState = dss.front();
235  const DetId idDetLay = detWithState.first->geographicalId();
236 
237  // check if this is a DT and the track has the needed quality
238  if(!chamberSelection(idDetLay,trans_track)) continue;
239 
240  DTChamberId DTid = (DTChamberId) idDetLay;
241 
242  // check if the chamber has already been counted
243  if(alreadyCheckedCh.find(DTid) != alreadyCheckedCh.end()) continue;
244  alreadyCheckedCh.insert(DTid);
245 
246  // get the compatible measurements
247  MeasurementContainer detMeasurements_initial = theMeasurementExtractor->measurements(layer_list.at(i),
248  detWithState.first,
249  detWithState.second,
250  *theEstimator, event);
251  LogTrace("DTDQM|DTMonitorModule|DTChamberEfficiency") << " chamber: " << DTid
252  << " has: " << detMeasurements_initial.size() << " comp. meas." << endl;
253 
254 
255  //we want to be more picky about the quality of the segments:
256  //exclude the segments with less than 12 hits
257  MeasurementContainer detMeasurements = segQualityCut(detMeasurements_initial);
258 
259  // get the histos for this chamber
260  vector<MonitorElement *> histos = histosPerW[DTid.wheel()+2];
261  // fill them
262  if (detMeasurements_initial.size() != 0) histos[0]->Fill(DTid.sector(),DTid.station(),1.);
263  if (detMeasurements.size() != 0) histos[1]->Fill(DTid.sector(),DTid.station(),1.);
264  histos[2]->Fill(DTid.sector(),DTid.station(),1.);
265 
266 
267  }
268  }
269  } else {
270  LogInfo("DTDQM|DTMonitorModule|DTChamberEfficiency") << "[DTChamberEfficiency] Collection: " << theTracksLabel
271  << " is not valid!" << endl;
272  }
273  return;
274 }
275 
276 bool DTChamberEfficiency::chamberSelection(const DetId& idDetLay, reco::TransientTrack& trans_track) const
277 {
278  //check that we have a muon and that is a DT detector
279  if(!(idDetLay.det() == DetId::Muon && idDetLay.subdetId() == MuonSubdetId::DT)) return false;
280 
281  if(trans_track.recHitsSize() == 2)
282  if(trans_track.recHit(0)->geographicalId() == idDetLay ||
283  trans_track.recHit(1)->geographicalId() == idDetLay) return false;
284 
285  return true;
286 }
287 
288 //riempi una per ogni segmento e una per segmento sopra 12 hit
289 
291 {
292 
294 
295  for(MeasurementContainer::const_iterator mescont_Itr = seg_list.begin();
296  mescont_Itr != seg_list.end(); ++mescont_Itr){
297 
298  //get the rechits of the segment
299  TransientTrackingRecHit::ConstRecHitContainer recHit_list = mescont_Itr->recHit()->transientHits();
300 
301  //loop over the rechits and get the number of hits
302  int nhit_seg(0);
303  for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator recList_Itr = recHit_list.begin();
304  recList_Itr != recHit_list.end(); ++recList_Itr){
305 
306  nhit_seg += (int)(*recList_Itr)->transientHits().size();
307  }
308 
309  DTChamberId tmpId = (DTChamberId)mescont_Itr->recHit()->hit()->geographicalId();
310 
311  if(tmpId.station() < 4 && nhit_seg >= 12) result.push_back(*mescont_Itr);
312  if(tmpId.station() == 4 && nhit_seg >= 8) result.push_back(*mescont_Itr);
313  }
314 
315  return result;
316 }
317 
318 vector<const DetLayer*> DTChamberEfficiency::compatibleLayers(const DetLayer *initialLayer,
319  const FreeTrajectoryState& fts, PropagationDirection propDir)
320 {
321 
322 vector<const DetLayer*> detLayers;
323 
324 if(theNavigationType == "Standard"){
325  // ask for compatible layers
326  detLayers = initialLayer->compatibleLayers(fts,propDir);
327  // I have to fit by hand the first layer until the seedTSOS is defined on the first rechit layer
328  // In fact the first layer is not returned by initialLayer->compatibleLayers.
329  detLayers.insert(detLayers.begin(),initialLayer);
330  }
331  else if (theNavigationType == "Direct"){
332 
333  DirectMuonNavigation navigation(ESHandle<MuonDetLayerGeometry>(&*theService->detLayerGeometry()));
334  detLayers = navigation.compatibleLayers(fts,propDir);
335  }
336  else
337  LogError("DTDQM|DTMonitorModule|DTChamberEfficiency") << "No Properly Navigation Selected!!"<<endl;
338 
339  return detLayers;
340 }
341 
343  return theService->propagator("SteppingHelixPropagatorAny");
344 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
DTChamberEfficiency(const edm::ParameterSet &pset)
size_t recHitsSize() const
number of RecHits
PropagationDirection
std::vector< const DetLayer * > compatibleLayers(const FreeTrajectoryState &fts, PropagationDirection timeDirection) const
void bookHistos()
Definition: Histogram.h:33
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
std::vector< const DetLayer * > compatibleLayers(const DetLayer *initialLayer, const FreeTrajectoryState &fts, PropagationDirection propDir)
bool chamberSelection(const DetId &idDetLay, reco::TransientTrack &trans_track) const
edm::ESHandle< Propagator > propagator() const
FreeTrajectoryState * freeState(bool withErrors=true) const
tuple result
Definition: query.py:137
TrackingRecHitRef recHit(size_t i) const
get n-th recHit
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
bool isValid() const
Definition: HandleBase.h:76
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
#define LogTrace(id)
void beginRun(const edm::Run &, const edm::EventSetup &)
std::vector< ConstRecHitPointer > ConstRecHitContainer
MeasurementContainer segQualityCut(const MeasurementContainer seg_list) const
Definition: DetId.h:20
dictionary histos
Definition: combine.py:3
tuple tracks
Definition: testEve_cfg.py:39
std::vector< TrajectoryMeasurement > MeasurementContainer
const T & get() const
Definition: EventSetup.h:55
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
int sector() const
Definition: DTChamberId.h:63
const Surface & surface() const
std::vector< const DetLayer * > compatibleLayers(NavigationDirection direction) const
Definition: DetLayer.cc:60
std::pair< const GeomDet *, TrajectoryStateOnSurface > DetWithState
static const int DT
Definition: MuonSubdetId.h:14
int station() const
Return the station number.
Definition: DTChamberId.h:53
#define debug
Definition: MEtoEDMFormat.h:34
Detector det() const
get the detector field from this detid
Definition: DetId.h:37
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:47
Definition: Run.h:31