CMS 3D CMS Logo

DTResidualCalibration.cc
Go to the documentation of this file.
1 
2 /*
3  * See header file for a description of this class.
4  *
5  */
6 
8 
9 // Framework
15 
16 //Geometry
19 
20 //RecHit
23 
27 
28 #include "TFile.h"
29 #include "TH1F.h"
30 #include "TH2F.h"
31 
32 #include <algorithm>
33 
35  histRange_(pset.getParameter<double>("histogramRange")),
36  segment4DLabel_(pset.getParameter<edm::InputTag>("segment4DLabel")),
37  rootBaseDir_(pset.getUntrackedParameter<std::string>("rootBaseDir","DT/Residuals")),
38  detailedAnalysis_(pset.getUntrackedParameter<bool>("detailedAnalysis",false)) {
39 
41  select_ = new DTSegmentSelector(pset,collector);
42 
43  LogDebug("Calibration") << "[DTResidualCalibration] Constructor called.";
44  consumes< DTRecSegment4DCollection >(edm::InputTag(segment4DLabel_));
45  std::string rootFileName = pset.getUntrackedParameter<std::string>("rootFileName","residuals.root");
46  rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
47  rootFile_->cd();
48 
49  segmok=0;
50  segmbad=0;
51  nevent=0;
52 }
53 
55  delete select_;
56  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Destructor called.";
57  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Analyzed events: " << nevent;
58  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Good segments: " << segmok;
59  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Bad segments: " << segmbad;
60 }
61 
63  TH1::SetDefaultSumw2(true);
64 }
65 
67 
68  // get the geometry
70  setup.get<MuonGeometryRecord>().get(dtGeomH);
71  dtGeom_ = dtGeomH.product();
72 
73  // Loop over all the chambers
74  if(histoMapTH1F_.empty()) {
75  for (auto ch_it : dtGeom_->chambers()) {
76  // Loop over the SLs
77  for (auto sl_it : ch_it->superLayers()) {
78  DTSuperLayerId slId = (sl_it)->id();
79  bookHistos(slId);
80  if(detailedAnalysis_) {
81  for (auto layer_it : (sl_it)->layers()) {
82  DTLayerId layerId = (layer_it)->id();
83  bookHistos(layerId);
84  }
85  }
86  }
87  }
88  }
89 }
90 
92 
93  rootFile_->cd();
94  ++nevent;
95 
96  // Get the 4D rechits from the event
98  event.getByLabel(segment4DLabel_, segments4D);
99 
100  // Loop over segments by chamber
102  for(chamberIdIt = segments4D->id_begin(); chamberIdIt != segments4D->id_end(); ++chamberIdIt){
103 
104  const DTChamber* chamber = dtGeom_->chamber(*chamberIdIt);
105 
106  // Get the range for the corresponding ChamberId
107  DTRecSegment4DCollection::range range = segments4D->get((*chamberIdIt));
108  // Loop over the rechits of this DetUnit
109  for(DTRecSegment4DCollection::const_iterator segment = range.first;
110  segment != range.second; ++segment){
111 
112  LogTrace("Calibration") << "Segment local pos (in chamber RF): " << (*segment).localPosition()
113  << "\nSegment global pos: " << chamber->toGlobal((*segment).localPosition());
114 
115  if( !(*select_)(*segment, event, setup) ) { segmbad++; continue; }
116  segmok++;
117 
118  // Get all 1D RecHits at step 3 within the 4D segment
119  std::vector<DTRecHit1D> recHits1D_S3;
120 
121  if( (*segment).hasPhi() ){
122  const DTChamberRecSegment2D* phiSeg = (*segment).phiSegment();
123  const std::vector<DTRecHit1D>& phiRecHits = phiSeg->specificRecHits();
124  std::copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
125  }
126 
127  if( (*segment).hasZed() ){
128  const DTSLRecSegment2D* zSeg = (*segment).zSegment();
129  const std::vector<DTRecHit1D>& zRecHits = zSeg->specificRecHits();
130  std::copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
131  }
132 
133  // Loop over 1D RecHit inside 4D segment
134  for(std::vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin();
135  recHit1D != recHits1D_S3.end(); ++recHit1D) {
136  const DTWireId wireId = recHit1D->wireId();
137 
138  float segmDistance = segmentToWireDistance(*recHit1D,*segment);
139  if(segmDistance > 2.1) LogTrace("Calibration") << "WARNING: segment-wire distance: " << segmDistance;
140  else LogTrace("Calibration") << "segment-wire distance: " << segmDistance;
141 
142  float residualOnDistance = DTRecHitSegmentResidual().compute(dtGeom_,*recHit1D,*segment);
143  LogTrace("Calibration") << "Wire Id " << wireId << " residual on distance: " << residualOnDistance;
144 
145  fillHistos(wireId.superlayerId(), segmDistance, residualOnDistance);
146  if(detailedAnalysis_) fillHistos(wireId.layerId(), segmDistance, residualOnDistance);
147  }
148  }
149  }
150 }
151 
153 
154  // Get the layer and the wire position
155  const DTWireId wireId = recHit1D.wireId();
156  const DTLayer* layer = dtGeom_->layer(wireId);
157  float wireX = layer->specificTopology().wirePosition(wireId.wire());
158 
159  // Extrapolate the segment to the z of the wire
160  // Get wire position in chamber RF
161  // (y and z must be those of the hit to be coherent in the transf. of RF in case of rotations of the layer alignment)
162  LocalPoint wirePosInLay(wireX,recHit1D.localPosition().y(),recHit1D.localPosition().z());
163  GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
164  const DTChamber* chamber = dtGeom_->chamber(wireId.layerId().chamberId());
165  LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
166 
167  // Segment position at Wire z in chamber local frame
168  LocalPoint segPosAtZWire = segment.localPosition() + segment.localDirection()*wirePosInChamber.z()/cos(segment.localDirection().theta());
169 
170  // Compute the distance of the segment from the wire
171  int sl = wireId.superlayer();
172  float segmDistance = -1;
173  if(sl == 1 || sl == 3) segmDistance = fabs(wirePosInChamber.x() - segPosAtZWire.x());
174  else if(sl == 2) segmDistance = fabs(segPosAtZWire.y() - wirePosInChamber.y());
175 
176  return segmDistance;
177 }
178 
180 
181  LogDebug("Calibration") << "[DTResidualCalibration] Writing histos to file.";
182  rootFile_->cd();
183  rootFile_->Write();
184  rootFile_->Close();
185 
186  /*std::map<DTSuperLayerId, TH1F* >::const_iterator itSlHistos = histoMapTH1F_.begin();
187  std::map<DTSuperLayerId, TH1F* >::const_iterator itSlHistos_end = histoMapTH1F_.end();
188  for(; itSlHistos != itSlHistos_end; ++itSlHistos){
189  std::vector<TH1F*>::const_iterator itHistTH1F = (*itSlHistos).second.begin();
190  std::vector<TH1F*>::const_iterator itHistTH1F_end = (*itSlHistos).second.end();
191  for(; itHistTH1F != itHistTH1F_end; ++itHistTH1F) (*itHistTH1F)->Write();
192 
193  std::vector<TH2F*>::const_iterator itHistTH2F = histoMapTH2F_[(*itSlHistos).first].begin();
194  std::vector<TH2F*>::const_iterator itHistTH2F_end = histoMapTH2F_[(*itSlHistos).first].end();
195  for(; itHistTH2F != itHistTH2F_end; ++itHistTH2F) (*itHistTH2F)->Write();
196  }*/
197 
198 }
199 
201  TH1AddDirectorySentry addDir;
202  rootFile_->cd();
203 
204  LogDebug("Calibration") << "[DTResidualCalibration] Booking histos for SL: " << slId;
205 
206  // Compose the chamber name
207  // Define the step
208  int step = 3;
209 
210  std::string wheelStr = std::to_string(slId.wheel());
211  std::string stationStr = std::to_string(slId.station());
212  std::string sectorStr = std::to_string(slId.sector());
213 
214  std::string slHistoName =
215  "_STEP" + std::to_string(step) +
216  "_W" + wheelStr +
217  "_St" + stationStr +
218  "_Sec" + sectorStr +
219  "_SL" + std::to_string(slId.superlayer());
220 
221  LogDebug("Calibration") << "Accessing " << rootBaseDir_;
222  TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
223  if(!baseDir) baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
224  LogDebug("Calibration") << "Accessing " << ("Wheel" + wheelStr);
225  TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr).c_str());
226  if(!wheelDir) wheelDir = baseDir->mkdir(("Wheel" + wheelStr).c_str());
227  LogDebug("Calibration") << "Accessing " << ("Station" + stationStr);
228  TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr).c_str());
229  if(!stationDir) stationDir = wheelDir->mkdir(("Station" + stationStr).c_str());
230  LogDebug("Calibration") << "Accessing " << ("Sector" + sectorStr);
231  TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr).c_str());
232  if(!sectorDir) sectorDir = stationDir->mkdir(("Sector" + sectorStr).c_str());
233 
234  sectorDir->cd();
235 
236  // Create the monitor elements
237  TH1F* histosTH1F = new TH1F(("hResDist"+slHistoName).c_str(),
238  "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
239  200, -histRange_, histRange_);
240  TH2F* histosTH2F = new TH2F(("hResDistVsDist"+slHistoName).c_str(),
241  "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
242  100, 0, 2.5, 200, -histRange_, histRange_);
243  histoMapTH1F_[slId] = histosTH1F;
244  histoMapTH2F_[slId] = histosTH2F;
245 }
246 
248  TH1AddDirectorySentry addDir;
249  rootFile_->cd();
250 
251  LogDebug("Calibration") << "[DTResidualCalibration] Booking histos for layer: " << layerId;
252 
253  // Compose the chamber name
254  std::string wheelStr = std::to_string(layerId.wheel());
255  std::string stationStr = std::to_string(layerId.station());
256  std::string sectorStr = std::to_string(layerId.sector());
257  std::string superLayerStr = std::to_string(layerId.superlayer());
258  std::string layerStr = std::to_string(layerId.layer());
259  // Define the step
260  int step = 3;
261 
262  std::string layerHistoName =
263  "_STEP" + std::to_string(step) +
264  "_W" + wheelStr +
265  "_St" + stationStr +
266  "_Sec" + sectorStr +
267  "_SL" + superLayerStr +
268  "_Layer" + layerStr;
269 
270  LogDebug("Calibration") << "Accessing " << rootBaseDir_;
271  TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
272  if(!baseDir) baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
273  LogDebug("Calibration") << "Accessing " << ("Wheel" + wheelStr);
274  TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr).c_str());
275  if(!wheelDir) wheelDir = baseDir->mkdir(("Wheel" + wheelStr).c_str());
276  LogDebug("Calibration") << "Accessing " << ("Station" + stationStr);
277  TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr).c_str());
278  if(!stationDir) stationDir = wheelDir->mkdir(("Station" + stationStr).c_str());
279  LogDebug("Calibration") << "Accessing " << ("Sector" + sectorStr);
280  TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr).c_str());
281  if(!sectorDir) sectorDir = stationDir->mkdir(("Sector" + sectorStr).c_str());
282  LogDebug("Calibration") << "Accessing " << ("SL" + superLayerStr);
283  TDirectory* superLayerDir = sectorDir->GetDirectory(("SL" + superLayerStr).c_str());
284  if(!superLayerDir) superLayerDir = sectorDir->mkdir(("SL" + superLayerStr).c_str());
285 
286  superLayerDir->cd();
287  // Create histograms
288  TH1F* histosTH1F = new TH1F(("hResDist"+layerHistoName).c_str(),
289  "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
290  200, -histRange_, histRange_);
291  TH2F* histosTH2F = new TH2F(("hResDistVsDist"+layerHistoName).c_str(),
292  "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
293  100, 0, 2.5, 200, -histRange_, histRange_);
294  histoMapPerLayerTH1F_[layerId] = histosTH1F;
295  histoMapPerLayerTH2F_[layerId] = histosTH2F;
296 }
297 
298 // Fill a set of histograms for a given SL
300  float distance,
301  float residualOnDistance) {
302  histoMapTH1F_[slId]->Fill(residualOnDistance);
303  histoMapTH2F_[slId]->Fill(distance, residualOnDistance);
304 }
305 
306 // Fill a set of histograms for a given layer
308  float distance,
309  float residualOnDistance) {
310  histoMapPerLayerTH1F_[layerId]->Fill(residualOnDistance);
311  histoMapPerLayerTH2F_[layerId]->Fill(distance, residualOnDistance);
312 }
313 
#define LogDebug(id)
std::map< DTSuperLayerId, TH1F * > histoMapTH1F_
LocalPoint localPosition() const override
Return the 3-dimensional local position.
Definition: DTRecHit1D.h:60
T getUntrackedParameter(std::string const &, T const &) const
float wirePosition(int wireNumber) const
Returns the x position in the layer of a given wire number.
Definition: DTTopology.cc:86
DTSegmentSelector * select_
const std::vector< const DTChamber * > & chambers() const
Return a vector of all Chamber.
Definition: DTGeometry.cc:102
LocalPoint localPosition() const override
Local position in Chamber frame.
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
Definition: LayerTriplets.cc:4
def copy(args, dbName)
void fillHistos(DTSuperLayerId slId, float distance, float residualOnDistance)
void bookHistos(DTSuperLayerId slId)
const DTChamber * chamber(const DTChamberId &id) const
Return a DTChamber given its id.
Definition: DTGeometry.cc:117
LocalVector localDirection() const override
Local direction in Chamber frame.
DTResidualCalibration(const edm::ParameterSet &pset)
Constructor.
DTChamberId chamberId() const
Return the corresponding ChamberId.
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
T y() const
Definition: PV3DBase.h:63
def setup(process, global_tag, zero_tesla=False)
Definition: GeneralSetup.py:2
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:69
identifier iterator
Definition: RangeMap.h:135
int layer() const
Return the layer number.
Definition: DTLayerId.h:53
DTSuperLayerId superlayerId() const
Return the corresponding SuperLayerId.
Definition: DTLayerId.h:59
std::map< DTLayerId, TH1F * > histoMapPerLayerTH1F_
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:43
const DTTopology & specificTopology() const
Definition: DTLayer.cc:42
void beginRun(const edm::Run &, const edm::EventSetup &) override
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::map< DTSuperLayerId, TH2F * > histoMapTH2F_
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
#define LogTrace(id)
std::vector< DTRecHit1D > specificRecHits() const
Access to specific components.
std::map< DTLayerId, TH2F * > histoMapPerLayerTH2F_
~DTResidualCalibration() override
Destructor.
int wire() const
Return the wire number.
Definition: DTWireId.h:56
int superlayer() const
Return the superlayer number (deprecated method name)
void analyze(const edm::Event &event, const edm::EventSetup &setup) override
float compute(const DTGeometry *, const DTRecHit1D &, const DTRecSegment4D &)
const DTGeometry * dtGeom_
DTLayerId layerId() const
Return the corresponding LayerId.
Definition: DTWireId.h:62
HLT enums.
int sector() const
Definition: DTChamberId.h:61
T get() const
Definition: EventSetup.h:71
step
Definition: StallMonitor.cc:94
const DTLayer * layer(const DTLayerId &id) const
Return a layer given its id.
Definition: DTGeometry.cc:127
float segmentToWireDistance(const DTRecHit1D &recHit1D, const DTRecSegment4D &segment)
int station() const
Return the station number.
Definition: DTChamberId.h:51
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:45
T x() const
Definition: PV3DBase.h:62
T const * product() const
Definition: ESHandle.h:86
Definition: event.py:1
Definition: Run.h:45
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:107