CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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  select_(pset),
36  segment4DLabel_(pset.getParameter<edm::InputTag>("segment4DLabel")),
37  rootBaseDir_(pset.getUntrackedParameter<std::string>("rootBaseDir","DT/Residuals")),
38  detailedAnalysis_(pset.getUntrackedParameter<bool>("detailedAnalysis",false)) {
39 
40  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Constructor called.";
41 
42  std::string rootFileName = pset.getUntrackedParameter<std::string>("rootFileName","residuals.root");
43  rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
44  rootFile_->cd();
45 }
46 
48  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Destructor called.";
49 }
50 
52  TH1::SetDefaultSumw2(true);
53 }
54 
56 
57  // get the geometry
59  setup.get<MuonGeometryRecord>().get(dtGeomH);
60  dtGeom_ = dtGeomH.product();
61 
62  // Loop over all the chambers
63  if(histoMapTH1F_.size() == 0) {
64  auto ch_it = dtGeom_->chambers().begin();
65  auto ch_end = dtGeom_->chambers().end();
66  for (; ch_it != ch_end; ++ch_it) {
67  std::vector<const DTSuperLayer*>::const_iterator sl_it = (*ch_it)->superLayers().begin();
68  std::vector<const DTSuperLayer*>::const_iterator sl_end = (*ch_it)->superLayers().end();
69  // Loop over the SLs
70  for(; sl_it != sl_end; ++sl_it) {
71  DTSuperLayerId slId = (*sl_it)->id();
72  bookHistos(slId);
73  if(detailedAnalysis_) {
74  std::vector<const DTLayer*>::const_iterator layer_it = (*sl_it)->layers().begin();
75  std::vector<const DTLayer*>::const_iterator layer_end = (*sl_it)->layers().end();
76  for(; layer_it != layer_end; ++layer_it) {
77  DTLayerId layerId = (*layer_it)->id();
78  bookHistos(layerId);
79  }
80  }
81  }
82  }
83  }
84 }
85 
87  rootFile_->cd();
88 
89  // Get the 4D rechits from the event
91  event.getByLabel(segment4DLabel_, segment4Ds);
92 
93  // Loop over segments by chamber
94  DTRecSegment4DCollection::id_iterator chamberIdIt;
95  for(chamberIdIt = segment4Ds->id_begin(); chamberIdIt != segment4Ds->id_end(); ++chamberIdIt){
96 
97  const DTChamber* chamber = dtGeom_->chamber(*chamberIdIt);
98 
99  // Get the range for the corresponding ChamberId
100  DTRecSegment4DCollection::range range = segment4Ds->get((*chamberIdIt));
101  // Loop over the rechits of this DetUnit
102  for(DTRecSegment4DCollection::const_iterator segment = range.first;
103  segment != range.second; ++segment){
104 
105  LogTrace("Calibration") << "Segment local pos (in chamber RF): " << (*segment).localPosition()
106  << "\nSegment global pos: " << chamber->toGlobal((*segment).localPosition());
107 
108  if( !select_(*segment, event, setup) ) continue;
109 
110  // Get all 1D RecHits at step 3 within the 4D segment
111  std::vector<DTRecHit1D> recHits1D_S3;
112 
113  if( (*segment).hasPhi() ){
114  const DTChamberRecSegment2D* phiSeg = (*segment).phiSegment();
115  const std::vector<DTRecHit1D>& phiRecHits = phiSeg->specificRecHits();
116  std::copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
117  }
118 
119  if( (*segment).hasZed() ){
120  const DTSLRecSegment2D* zSeg = (*segment).zSegment();
121  const std::vector<DTRecHit1D>& zRecHits = zSeg->specificRecHits();
122  std::copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
123  }
124 
125  // Loop over 1D RecHit inside 4D segment
126  for(std::vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin();
127  recHit1D != recHits1D_S3.end(); ++recHit1D) {
128  const DTWireId wireId = recHit1D->wireId();
129 
130  float segmDistance = segmentToWireDistance(*recHit1D,*segment);
131  if(segmDistance > 2.1) LogTrace("Calibration") << "WARNING: segment-wire distance: " << segmDistance;
132  else LogTrace("Calibration") << "segment-wire distance: " << segmDistance;
133 
134  float residualOnDistance = DTRecHitSegmentResidual().compute(dtGeom_,*recHit1D,*segment);
135  LogTrace("Calibration") << "Wire Id " << wireId << " residual on distance: " << residualOnDistance;
136 
137  fillHistos(wireId.superlayerId(), segmDistance, residualOnDistance);
138  if(detailedAnalysis_) fillHistos(wireId.layerId(), segmDistance, residualOnDistance);
139  }
140  }
141  }
142 
143 }
144 
146 
147  // Get the layer and the wire position
148  const DTWireId wireId = recHit1D.wireId();
149  const DTLayer* layer = dtGeom_->layer(wireId);
150  float wireX = layer->specificTopology().wirePosition(wireId.wire());
151 
152  // Extrapolate the segment to the z of the wire
153  // Get wire position in chamber RF
154  // (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)
155  LocalPoint wirePosInLay(wireX,recHit1D.localPosition().y(),recHit1D.localPosition().z());
156  GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
157  const DTChamber* chamber = dtGeom_->chamber(wireId.layerId().chamberId());
158  LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
159 
160  // Segment position at Wire z in chamber local frame
161  LocalPoint segPosAtZWire = segment.localPosition() + segment.localDirection()*wirePosInChamber.z()/cos(segment.localDirection().theta());
162 
163  // Compute the distance of the segment from the wire
164  int sl = wireId.superlayer();
165  float segmDistance = -1;
166  if(sl == 1 || sl == 3) segmDistance = fabs(wirePosInChamber.x() - segPosAtZWire.x());
167  else if(sl == 2) segmDistance = fabs(segPosAtZWire.y() - wirePosInChamber.y());
168 
169  return segmDistance;
170 }
171 
173 
174  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Writing histos to file.";
175  rootFile_->cd();
176  rootFile_->Write();
177  rootFile_->Close();
178 
179  /*std::map<DTSuperLayerId, std::vector<TH1F*> >::const_iterator itSlHistos = histoMapTH1F_.begin();
180  std::map<DTSuperLayerId, std::vector<TH1F*> >::const_iterator itSlHistos_end = histoMapTH1F_.end();
181  for(; itSlHistos != itSlHistos_end; ++itSlHistos){
182  std::vector<TH1F*>::const_iterator itHistTH1F = (*itSlHistos).second.begin();
183  std::vector<TH1F*>::const_iterator itHistTH1F_end = (*itSlHistos).second.end();
184  for(; itHistTH1F != itHistTH1F_end; ++itHistTH1F) (*itHistTH1F)->Write();
185 
186  std::vector<TH2F*>::const_iterator itHistTH2F = histoMapTH2F_[(*itSlHistos).first].begin();
187  std::vector<TH2F*>::const_iterator itHistTH2F_end = histoMapTH2F_[(*itSlHistos).first].end();
188  for(; itHistTH2F != itHistTH2F_end; ++itHistTH2F) (*itHistTH2F)->Write();
189  }*/
190 
191 }
192 
194  TH1AddDirectorySentry addDir;
195  rootFile_->cd();
196 
197  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Booking histos for SL: " << slId;
198 
199  // Compose the chamber name
200  std::stringstream wheelStr; wheelStr << slId.wheel();
201  std::stringstream stationStr; stationStr << slId.station();
202  std::stringstream sectorStr; sectorStr << slId.sector();
203  std::stringstream superLayerStr; superLayerStr << slId.superlayer();
204  // Define the step
205  int step = 3;
206  std::stringstream stepStr; stepStr << step;
207 
208  std::string slHistoName =
209  "_STEP" + stepStr.str() +
210  "_W" + wheelStr.str() +
211  "_St" + stationStr.str() +
212  "_Sec" + sectorStr.str() +
213  "_SL" + superLayerStr.str();
214 
215  edm::LogVerbatim("Calibration") << "Accessing " << rootBaseDir_;
216  TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
217  if(!baseDir) baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
218  edm::LogVerbatim("Calibration") << "Accessing " << ("Wheel" + wheelStr.str());
219  TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr.str()).c_str());
220  if(!wheelDir) wheelDir = baseDir->mkdir(("Wheel" + wheelStr.str()).c_str());
221  edm::LogVerbatim("Calibration") << "Accessing " << ("Station" + stationStr.str());
222  TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr.str()).c_str());
223  if(!stationDir) stationDir = wheelDir->mkdir(("Station" + stationStr.str()).c_str());
224  edm::LogVerbatim("Calibration") << "Accessing " << ("Sector" + sectorStr.str());
225  TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr.str()).c_str());
226  if(!sectorDir) sectorDir = stationDir->mkdir(("Sector" + sectorStr.str()).c_str());
227 
228  /*std::string dirName = rootBaseDir_ + "/Wheel" + wheelStr.str() +
229  "/Station" + stationStr.str() +
230  "/Sector" + sectorStr.str();
231 
232  TDirectory* dir = rootFile_->GetDirectory(dirName.c_str());
233  if(!dir) dir = rootFile_->mkdir(dirName.c_str());
234  dir->cd();*/
235  sectorDir->cd();
236  // Create the monitor elements
237  std::vector<TH1F*> histosTH1F;
238  histosTH1F.push_back(new TH1F(("hResDist"+slHistoName).c_str(),
239  "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
240  200, -0.4, 0.4));
241  std::vector<TH2F*> histosTH2F;
242  histosTH2F.push_back(new TH2F(("hResDistVsDist"+slHistoName).c_str(),
243  "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
244  100, 0, 2.5, 200, -0.4, 0.4));
245  histoMapTH1F_[slId] = histosTH1F;
246  histoMapTH2F_[slId] = histosTH2F;
247 }
248 
250  TH1AddDirectorySentry addDir;
251  rootFile_->cd();
252 
253  edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Booking histos for layer: " << layerId;
254 
255  // Compose the chamber name
256  std::stringstream wheelStr; wheelStr << layerId.wheel();
257  std::stringstream stationStr; stationStr << layerId.station();
258  std::stringstream sectorStr; sectorStr << layerId.sector();
259  std::stringstream superLayerStr; superLayerStr << layerId.superlayer();
260  std::stringstream layerStr; layerStr << layerId.layer();
261  // Define the step
262  int step = 3;
263  std::stringstream stepStr; stepStr << step;
264 
265  std::string layerHistoName =
266  "_STEP" + stepStr.str() +
267  "_W" + wheelStr.str() +
268  "_St" + stationStr.str() +
269  "_Sec" + sectorStr.str() +
270  "_SL" + superLayerStr.str() +
271  "_Layer" + layerStr.str();
272 
273  edm::LogVerbatim("Calibration") << "Accessing " << rootBaseDir_;
274  TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
275  if(!baseDir) baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
276  edm::LogVerbatim("Calibration") << "Accessing " << ("Wheel" + wheelStr.str());
277  TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr.str()).c_str());
278  if(!wheelDir) wheelDir = baseDir->mkdir(("Wheel" + wheelStr.str()).c_str());
279  edm::LogVerbatim("Calibration") << "Accessing " << ("Station" + stationStr.str());
280  TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr.str()).c_str());
281  if(!stationDir) stationDir = wheelDir->mkdir(("Station" + stationStr.str()).c_str());
282  edm::LogVerbatim("Calibration") << "Accessing " << ("Sector" + sectorStr.str());
283  TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr.str()).c_str());
284  if(!sectorDir) sectorDir = stationDir->mkdir(("Sector" + sectorStr.str()).c_str());
285  edm::LogVerbatim("Calibration") << "Accessing " << ("SL" + superLayerStr.str());
286  TDirectory* superLayerDir = sectorDir->GetDirectory(("SL" + superLayerStr.str()).c_str());
287  if(!superLayerDir) superLayerDir = sectorDir->mkdir(("SL" + superLayerStr.str()).c_str());
288 
289  superLayerDir->cd();
290  // Create histograms
291  std::vector<TH1F*> histosTH1F;
292  histosTH1F.push_back(new TH1F(("hResDist"+layerHistoName).c_str(),
293  "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
294  200, -0.4, 0.4));
295  std::vector<TH2F*> histosTH2F;
296  histosTH2F.push_back(new TH2F(("hResDistVsDist"+layerHistoName).c_str(),
297  "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
298  100, 0, 2.5, 200, -0.4, 0.4));
299  histoMapPerLayerTH1F_[layerId] = histosTH1F;
300  histoMapPerLayerTH2F_[layerId] = histosTH2F;
301 }
302 
303 // Fill a set of histograms for a given SL
305  float distance,
306  float residualOnDistance) {
307  std::vector<TH1F*> const& histosTH1F = histoMapTH1F_[slId];
308  std::vector<TH2F*> const& histosTH2F = histoMapTH2F_[slId];
309  histosTH1F[0]->Fill(residualOnDistance);
310  histosTH2F[0]->Fill(distance, residualOnDistance);
311 }
312 
313 // Fill a set of histograms for a given layer
315  float distance,
316  float residualOnDistance) {
317  std::vector<TH1F*> const& histosTH1F = histoMapPerLayerTH1F_[layerId];
318  std::vector<TH2F*> const& histosTH2F = histoMapPerLayerTH2F_[layerId];
319  histosTH1F[0]->Fill(residualOnDistance);
320  histosTH2F[0]->Fill(distance, residualOnDistance);
321 }
322 
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
const std::vector< const DTChamber * > & chambers() const
Return a vector of all Chamber.
Definition: DTGeometry.cc:85
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
void fillHistos(DTSuperLayerId slId, float distance, float residualOnDistance)
void bookHistos(DTSuperLayerId slId)
std::map< DTLayerId, std::vector< TH1F * > > histoMapPerLayerTH1F_
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
virtual ~DTResidualCalibration()
Destructor.
T y() const
Definition: PV3DBase.h:63
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:69
int layer() const
Return the layer number.
Definition: DTLayerId.h:53
DTSuperLayerId superlayerId() const
Return the corresponding SuperLayerId.
Definition: DTLayerId.h:59
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
void analyze(const edm::Event &event, const edm::EventSetup &setup)
const DTLayer * layer(DTLayerId id) const
Return a layer given its id.
Definition: DTGeometry.cc:110
virtual LocalVector localDirection() const
Local direction in Chamber frame.
void beginRun(const edm::Run &, const edm::EventSetup &)
const DTTopology & specificTopology() const
Definition: DTLayer.cc:42
std::map< DTSuperLayerId, std::vector< TH1F * > > histoMapTH1F_
std::map< DTLayerId, std::vector< TH2F * > > histoMapPerLayerTH2F_
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::map< DTSuperLayerId, std::vector< 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
virtual LocalPoint localPosition() const
Local position in Chamber frame.
#define LogTrace(id)
std::vector< DTRecHit1D > specificRecHits() const
Access to specific components.
virtual LocalPoint localPosition() const
Return the 3-dimensional local position.
Definition: DTRecHit1D.h:60
int wire() const
Return the wire number.
Definition: DTWireId.h:56
int superlayer() const
Return the superlayer number (deprecated method name)
float compute(const DTGeometry *, const DTRecHit1D &, const DTRecSegment4D &)
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
const DTChamber * chamber(DTChamberId id) const
Return a DTChamber given its id.
Definition: DTGeometry.cc:100
const DTGeometry * dtGeom_
DTLayerId layerId() const
Return the corresponding LayerId.
Definition: DTWireId.h:62
int sector() const
Definition: DTChamberId.h:61
volatile std::atomic< bool > shutdown_flag false
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
Definition: Run.h:43
DTWireId wireId() const
Return the wireId.
Definition: DTRecHit1D.h:107