CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DTSegment4DQuality.cc
Go to the documentation of this file.
1 /*
2  * See header file for a description of this class.
3  *
4  * $Date: 2009/11/04 17:22:32 $
5  * $Revision: 1.9 $
6  * \author S. Bolognesi and G. Cerminara - INFN Torino
7  */
8 
9 #include "DTSegment4DQuality.h"
11 
16 
20 
27 
28 #include "Histograms.h"
29 
30 #include "TFile.h"
31 #include <iostream>
32 #include <map>
33 
34 
35 using namespace std;
36 using namespace edm;
37 
38 // Book the histos
39 
40 // Constructor
42  // Get the debug parameter for verbose output
43  debug = pset.getUntrackedParameter<bool>("debug");
45 
46  rootFileName = pset.getUntrackedParameter<string>("rootFileName");
47  // the name of the simhit collection
48  simHitLabel = pset.getUntrackedParameter<InputTag>("simHitLabel");
49  // the name of the 4D rec hit collection
50  segment4DLabel = pset.getUntrackedParameter<InputTag>("segment4DLabel");
51 
52  //sigma resolution on position
53  sigmaResX = pset.getParameter<double>("sigmaResX");
54  sigmaResY = pset.getParameter<double>("sigmaResY");
55  //sigma resolution on angle
56  sigmaResAlpha = pset.getParameter<double>("sigmaResAlpha");
57  sigmaResBeta = pset.getParameter<double>("sigmaResBeta");
58  doall = pset.getUntrackedParameter<bool>("doall", false);
59  local = pset.getUntrackedParameter<bool>("local", false);
60 
61 
62  // Create the root file
63  //theFile = new TFile(rootFileName.c_str(), "RECREATE");
64  //theFile->cd();
65 // ----------------------
66  // get hold of back-end interface
67  dbe_ = 0;
68  dbe_ = Service<DQMStore>().operator->();
69  if ( dbe_ ) {
70  if (debug) {
71  dbe_->setVerbose(1);
72  } else {
73  dbe_->setVerbose(0);
74  }
75  }
76  if ( dbe_ ) {
77  if ( debug ) dbe_->showDirStructure();
78  }
79 
80  h4DHit= new HRes4DHit ("All",dbe_,doall,local);
81  h4DHit_W0= new HRes4DHit ("W0",dbe_,doall,local);
82  h4DHit_W1= new HRes4DHit ("W1",dbe_,doall,local);
83  h4DHit_W2= new HRes4DHit ("W2",dbe_,doall,local);
84 
85  if(doall) {
86  hEff_All= new HEff4DHit ("All",dbe_);
87  hEff_W0= new HEff4DHit ("W0",dbe_);
88  hEff_W1= new HEff4DHit ("W1",dbe_);
89  hEff_W2= new HEff4DHit ("W2",dbe_);
90  }
91 }
92 
93 // Destructor
95 
96 }
98  edm::EventSetup const& c){
99 
100 
101 }
102 
104  // Write the histos to file
105  //theFile->cd();
106 
107  //h4DHit->Write();
108  //h4DHit_W0->Write();
109  //h4DHit_W1->Write();
110  //h4DHit_W2->Write();
111 
112  if(doall){
113  hEff_All->ComputeEfficiency();
114  hEff_W0->ComputeEfficiency();
115  hEff_W1->ComputeEfficiency();
116  hEff_W2->ComputeEfficiency();
117  }
118  //hEff_All->Write();
119  //hEff_W0->Write();
120  //hEff_W1->Write();
121  //hEff_W2->Write();
122  //if ( rootFileName.size() != 0 && dbe_ ) dbe_->save(rootFileName);
123 
124  //theFile->Close();
125 }
126 
127 // The real analysis
128  void DTSegment4DQuality::analyze(const Event & event, const EventSetup& eventSetup){
129  //theFile->cd();
130 
131  // Get the DT Geometry
132  ESHandle<DTGeometry> dtGeom;
133  eventSetup.get<MuonGeometryRecord>().get(dtGeom);
134 
135  // Get the SimHit collection from the event
137  event.getByLabel(simHitLabel, simHits); //FIXME: second string to be removed
138 
139  //Map simHits by chamber
140  map<DTChamberId, PSimHitContainer > simHitsPerCh;
141  for(PSimHitContainer::const_iterator simHit = simHits->begin();
142  simHit != simHits->end(); simHit++){
143  // Create the id of the chamber (the simHits in the DT known their wireId)
144  DTChamberId chamberId = (((DTWireId(simHit->detUnitId())).layerId()).superlayerId()).chamberId();
145  // Fill the map
146  simHitsPerCh[chamberId].push_back(*simHit);
147  }
148 
149  // Get the 4D rechits from the event
151  event.getByLabel(segment4DLabel, segment4Ds);
152 
153  if(!segment4Ds.isValid()) {
154  if(debug) cout << "[DTSegment4DQuality]**Warning: no 4D Segments with label: " <<segment4DLabel
155  << " in this event, skipping!" << endl;
156  return;
157  }
158 
159  // Loop over all chambers containing a segment
161  for (chamberId = segment4Ds->id_begin();
162  chamberId != segment4Ds->id_end();
163  ++chamberId){
164 
165  if((*chamberId).station() == 4)
166  continue; //use DTSegment2DQuality to analyze MB4 performaces
167 
168  //------------------------- simHits ---------------------------//
169  //Get simHits of each chamber
170  PSimHitContainer simHits = simHitsPerCh[(*chamberId)];
171 
172  // Map simhits per wire
173  map<DTWireId, PSimHitContainer > simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
174  map<DTWireId, const PSimHit*> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
175  int nMuSimHit = muSimHitPerWire.size();
176  if(nMuSimHit == 0 || nMuSimHit == 1) {
177  if(debug && nMuSimHit == 1)
178  cout << "[DTSegment4DQuality] Only " << nMuSimHit << " mu SimHit in this chamber, skipping!" << endl;
179  continue; // If no or only one mu SimHit is found skip this chamber
180  }
181  if(debug)
182  cout << "=== Chamber " << (*chamberId) << " has " << nMuSimHit << " SimHits" << endl;
183 
184  //Find outer and inner mu SimHit to build a segment
185  pair<const PSimHit*, const PSimHit*> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
186 
187  //Find direction and position of the sim Segment in Chamber RF
188  pair<LocalVector, LocalPoint> dirAndPosSimSegm = DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit,
189  (*chamberId),&(*dtGeom));
190 
191  LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
192  LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
193  const DTChamber* chamber = dtGeom->chamber(*chamberId);
194  GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos);
195  GlobalVector simSegmGlobalDir = chamber->toGlobal(simSegmLocalDir);
196 
197  //phi and theta angle of simulated segment in Chamber RF
198  float alphaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
199  float betaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).second;
200  //x,y position of simulated segment in Chamber RF
201  float xSimSeg = simSegmLocalPos.x();
202  float ySimSeg = simSegmLocalPos.y();
203  //Position (in eta,phi coordinates) in lobal RF
204  float etaSimSeg = simSegmGlobalPos.eta();
205  float phiSimSeg = simSegmGlobalPos.phi();
206 
207  if(debug)
208  cout<<" Simulated segment: local direction "<<simSegmLocalDir<<endl
209  <<" local position "<<simSegmLocalPos<<endl
210  <<" alpha "<<alphaSimSeg<<endl
211  <<" beta "<<betaSimSeg<<endl;
212 
213  //---------------------------- recHits --------------------------//
214  // Get the range of rechit for the corresponding chamberId
215  bool recHitFound = false;
216  DTRecSegment4DCollection::range range = segment4Ds->get(*chamberId);
217  int nsegm = distance(range.first, range.second);
218  if(debug)
219  cout << " Chamber: " << *chamberId << " has " << nsegm
220  << " 4D segments" << endl;
221 
222  if (nsegm!=0) {
223  // Find the best RecHit: look for the 4D RecHit with the phi angle closest
224  // to that of segment made of SimHits.
225  // RecHits must have delta alpha and delta position within 5 sigma of
226  // the residual distribution (we are looking for residuals of segments
227  // usefull to the track fit) for efficency purpose
228  const DTRecSegment4D* bestRecHit = 0;
229  bool bestRecHitFound = false;
230  double deltaAlpha = 99999;
231 
232  // Loop over the recHits of this chamberId
233  for (DTRecSegment4DCollection::const_iterator segment4D = range.first;
234  segment4D!=range.second;
235  ++segment4D){
236  // Check the dimension
237  if((*segment4D).dimension() != 4) {
238  if(debug)cout << "[DTSegment4DQuality]***Error: This is not 4D segment!!!" << endl;
239  continue;
240  }
241  // Segment Local Direction and position (in Chamber RF)
242  LocalVector recSegDirection = (*segment4D).localDirection();
243  //LocalPoint recSegPosition = (*segment4D).localPosition();
244 
245  float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
246  if(debug)
247  cout << " RecSegment direction: " << recSegDirection << endl
248  << " position : " << (*segment4D).localPosition() << endl
249  << " alpha : " << recSegAlpha << endl
250  << " beta : " << DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).second << endl;
251 
252  if(fabs(recSegAlpha - alphaSimSeg) < deltaAlpha) {
253  deltaAlpha = fabs(recSegAlpha - alphaSimSeg);
254  bestRecHit = &(*segment4D);
255  bestRecHitFound = true;
256  }
257  } // End of Loop over all 4D RecHits
258 
259  if(bestRecHitFound) {
260  // Best rechit direction and position in Chamber RF
261  LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
262  LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
263  // Errors on x and y
264  LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
265  LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
266 
267  float alphaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
268  float betaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).second;
269  // Errors on alpha and beta
270 
271  // Get position and direction using the rx projection (so in SL
272  // reference frame). Note that x (and y) are swapped wrt to Chamber
273  // frame
274  //if (bestRecHit->hasZed() ) {
275  const DTSLRecSegment2D * zedRecSeg = bestRecHit->zSegment();
276  LocalPoint bestRecHitLocalPosRZ = zedRecSeg->localPosition();
277  LocalVector bestRecHitLocalDirRZ = zedRecSeg->localDirection();
278  // Errors on x and y
279  LocalError bestRecHitLocalPosErrRZ = zedRecSeg->localPositionError();
280  LocalError bestRecHitLocalDirErrRZ = zedRecSeg->localDirectionError();
281 
282  float alphaBestRHitRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDirRZ).first;
283 
284  // Get SimSeg position and Direction in rZ SL frame
285  const DTSuperLayer* sl = dtGeom->superLayer(zedRecSeg->superLayerId());
286  LocalPoint simSegLocalPosRZTmp = sl->toLocal(simSegmGlobalPos);
287  LocalVector simSegLocalDirRZ = sl->toLocal(simSegmGlobalDir);
288  LocalPoint simSegLocalPosRZ =
289  simSegLocalPosRZTmp + simSegLocalDirRZ*(-simSegLocalPosRZTmp.z()/(cos(simSegLocalDirRZ.theta())));
290  float alphaSimSegRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegLocalDirRZ).first;
291 
292  if (debug) cout <<
293  "RZ SL: recPos " << bestRecHitLocalPosRZ <<
294  "recDir " << bestRecHitLocalDirRZ <<
295  "recAlpha " << alphaBestRHitRZ << endl <<
296  "RZ SL: simPos " << simSegLocalPosRZ <<
297  "simDir " << simSegLocalDirRZ <<
298  "simAlpha " << alphaSimSegRZ << endl ;
299  //}
300 
301 
302  if(fabs(alphaBestRHit - alphaSimSeg) < 5*sigmaResAlpha &&
303  fabs(betaBestRHit - betaSimSeg) < 5*sigmaResBeta &&
304  fabs(bestRecHitLocalPos.x() - xSimSeg) < 5*sigmaResX &&
305  fabs(bestRecHitLocalPos.y() - ySimSeg) < 5*sigmaResY) {
306  recHitFound = true;
307  }
308 
309  // Fill Residual histos
310  HRes4DHit *histo=0;
311 
312  if((*chamberId).wheel() == 0)
313  histo = h4DHit_W0;
314  else if(abs((*chamberId).wheel()) == 1)
315  histo = h4DHit_W1;
316  else if(abs((*chamberId).wheel()) == 2)
317  histo = h4DHit_W2;
318 
319  histo->Fill(alphaSimSeg,
320  alphaBestRHit,
321  betaSimSeg,
322  betaBestRHit,
323  xSimSeg,
324  bestRecHitLocalPos.x(),
325  ySimSeg,
326  bestRecHitLocalPos.y(),
327  etaSimSeg,
328  phiSimSeg,
329  bestRecHitLocalPosRZ.x(),
330  simSegLocalPosRZ.x(),
331  alphaBestRHitRZ,
332  alphaSimSegRZ,
333  sqrt(bestRecHitLocalDirErr.xx()),
334  sqrt(bestRecHitLocalDirErr.yy()),
335  sqrt(bestRecHitLocalPosErr.xx()),
336  sqrt(bestRecHitLocalPosErr.yy()),
337  sqrt(bestRecHitLocalDirErrRZ.xx()),
338  sqrt(bestRecHitLocalPosErrRZ.xx())
339  );
340 
341  h4DHit->Fill(alphaSimSeg,
342  alphaBestRHit,
343  betaSimSeg,
344  betaBestRHit,
345  xSimSeg,
346  bestRecHitLocalPos.x(),
347  ySimSeg,
348  bestRecHitLocalPos.y(),
349  etaSimSeg,
350  phiSimSeg,
351  bestRecHitLocalPosRZ.x(),
352  simSegLocalPosRZ.x(),
353  alphaBestRHitRZ,
354  alphaSimSegRZ,
355  sqrt(bestRecHitLocalDirErr.xx()),
356  sqrt(bestRecHitLocalDirErr.yy()),
357  sqrt(bestRecHitLocalPosErr.xx()),
358  sqrt(bestRecHitLocalPosErr.yy()),
359  sqrt(bestRecHitLocalDirErrRZ.xx()),
360  sqrt(bestRecHitLocalPosErrRZ.xx())
361  );
362  }
363  } //end of if(nsegm!=0)
364 
365  // Fill Efficiency plot
366  if(doall){
367  HEff4DHit *heff = 0;
368 
369  if((*chamberId).wheel() == 0)
370  heff = hEff_W0;
371  else if(abs((*chamberId).wheel()) == 1)
372  heff = hEff_W1;
373  else if(abs((*chamberId).wheel()) == 2)
374  heff = hEff_W2;
375  heff->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound);
376  hEff_All->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound);
377  }
378  } // End of loop over chambers
379  }
T getParameter(std::string const &) const
virtual LocalError localPositionError() const
local position error in SL frame
T getUntrackedParameter(std::string const &, T const &) const
float xx() const
Definition: LocalError.h:24
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:52
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:47
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
#define abs(x)
Definition: mlp_lapack.h:159
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:62
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
Perform the real analysis.
identifier iterator
Definition: RangeMap.h:138
virtual LocalError localDirectionError() const
the local direction error (xx,xy,yy) in SL frame: only xx is not 0.
virtual ~DTSegment4DQuality()
Destructor.
static std::map< DTWireId, const PSimHit * > mapMuSimHitsPerWire(const std::map< DTWireId, edm::PSimHitContainer > &simHitWireMap)
Create a map between the Mu SimHits and corresponding MuBarWireId ;.
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
virtual LocalError localPositionError() const
Local position error in Chamber frame.
static std::pair< const PSimHit *, const PSimHit * > findMuSimSegment(const std::map< DTWireId, const PSimHit * > &mapWireAndMuSimHit)
Find Innermost and outermost SimHit from Mu in a SL (they identify a simulated segment) ...
virtual LocalVector localDirection() const
Local direction in Chamber frame.
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:45
DTSegment4DQuality(const edm::ParameterSet &pset)
Constructor.
float yy() const
Definition: LocalError.h:26
T sqrt(T t)
Definition: SSEVec.h:48
void Fill(float simDirectionAlpha, float recDirectionAlpha, float simDirectionBeta, float recDirectionBeta, float simX, float recX, float simY, float recY, float simEta, float simPhi, float recYRZ, float simYRZ, float recBetaRZ, float simBetaRZ, float sigmaAlpha, float sigmaBeta, float sigmaX, float sigmaY, float sigmaBetaRZ, float sigmaYRZ)
Definition: Histograms.h:907
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
A set of histograms for efficiency 4D RecHits.
Definition: Histograms.h:1120
void endLuminosityBlock(edm::LuminosityBlock const &lumiSeg, edm::EventSetup const &c)
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
virtual LocalPoint localPosition() const
Local position in Chamber frame.
static std::pair< double, double > findSegmentAlphaAndBeta(const LocalVector &direction)
Find the angles from a segment direction:
DQMStore * dbe_
DTSuperLayerId superLayerId() const
The id of the superlayer on which reside the segment.
virtual LocalPoint localPosition() const
local position in SL frame
const DTSLRecSegment2D * zSegment() const
The Z segment: 0 if not zed projection available.
static std::pair< LocalVector, LocalPoint > findMuSimSegmentDirAndPos(const std::pair< const PSimHit *, const PSimHit * > &inAndOutSimHit, const DetId detId, const DTGeometry *muonGeom)
Find direction and position of a segment (in local RF) from outer and inner mu SimHit in the RF of ob...
const T & get() const
Definition: EventSetup.h:55
tuple simHits
Definition: trackerHits.py:16
T eta() const
Definition: PV3DBase.h:76
virtual LocalVector localDirection() const
the local direction in SL frame
tuple cout
Definition: gather_cfg.py:121
std::vector< PSimHit > PSimHitContainer
static std::map< DTWireId, edm::PSimHitContainer > mapSimHitsPerWire(const edm::PSimHitContainer &simhits)
#define debug
Definition: MEtoEDMFormat.h:34
T x() const
Definition: PV3DBase.h:62
void Fill(float etaSimSegm, float phiSimSegm, float xSimSegm, float ySimSegm, float alphaSimSegm, float betaSimSegm, bool fillRecHit)
Definition: Histograms.h:1215
virtual LocalError localDirectionError() const
Local direction error in the Chamber frame.