CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonAlignment.cc
Go to the documentation of this file.
1 
2 /*
3  * DQM muon alignment analysis monitoring
4  *
5  * \author J. Fernandez - Univ. Oviedo <Javier.Fernandez@cern.ch>
6  */
7 
9 
10 
12 
13  metname = "MuonAlignment";
14 
15  LogTrace(metname)<<"[MuonAlignment] Constructor called!"<<std::endl;
16 
17  parameters = pSet;
18 
19  theMuonCollectionLabel = consumes<reco::TrackCollection>(parameters.getParameter<edm::InputTag>("MuonCollection"));
20 
21  theRecHits4DTagDT = consumes<DTRecSegment4DCollection>(parameters.getParameter<edm::InputTag>("RecHits4DDTCollectionTag"));
22  theRecHits4DTagCSC = consumes<CSCSegmentCollection>(parameters.getParameter<edm::InputTag>("RecHits4DCSCCollectionTag"));
23 
24  resLocalXRangeStation1 = parameters.getUntrackedParameter<double>("resLocalXRangeStation1");
25  resLocalXRangeStation2 = parameters.getUntrackedParameter<double>("resLocalXRangeStation2");
26  resLocalXRangeStation3 = parameters.getUntrackedParameter<double>("resLocalXRangeStation3");
27  resLocalXRangeStation4 = parameters.getUntrackedParameter<double>("resLocalXRangeStation4");
28  resLocalYRangeStation1 = parameters.getUntrackedParameter<double>("resLocalYRangeStation1");
29  resLocalYRangeStation2 = parameters.getUntrackedParameter<double>("resLocalYRangeStation2");
30  resLocalYRangeStation3 = parameters.getUntrackedParameter<double>("resLocalYRangeStation3");
31  resLocalYRangeStation4 = parameters.getUntrackedParameter<double>("resLocalYRangeStation4");
32  resPhiRange = parameters.getUntrackedParameter<double>("resPhiRange");
33  resThetaRange = parameters.getUntrackedParameter<double>("resThetaRange");
34 
35  meanPositionRange = parameters.getUntrackedParameter<double>("meanPositionRange");
36  rmsPositionRange = parameters.getUntrackedParameter<double>("rmsPositionRange");
37  meanAngleRange = parameters.getUntrackedParameter<double>("meanAngleRange");
38  rmsAngleRange = parameters.getUntrackedParameter<double>("rmsAngleRange");
39 
40  nbins = parameters.getUntrackedParameter<unsigned int>("nbins");
41  min1DTrackRecHitSize = parameters.getUntrackedParameter<unsigned int>("min1DTrackRecHitSize");
42  min4DTrackSegmentSize = parameters.getUntrackedParameter<unsigned int>("min4DTrackSegmentSize");
43 
44  doDT = parameters.getUntrackedParameter<bool>("doDT");
45  doCSC = parameters.getUntrackedParameter<bool>("doCSC");
46  doSummary = parameters.getUntrackedParameter<bool>("doSummary");
47 
49  numberOfHits=0;
50 
52  topFolder << MEFolderName+"/Alignment/Muon";
53 }
54 
56 }
57 
58 
60 
61 
62  LogTrace(metname)<<"[MuonAlignment] Parameters initialization";
63 
64  if(!(doDT || doCSC) ) {
65  edm::LogError("MuonAlignment") <<" Error!! At least one Muon subsystem (DT or CSC) must be monitorized!!" << std::endl;
66  edm::LogError("MuonAlignment") <<" Please enable doDT or doCSC to True in your python cfg file!!!" << std::endl;
67  exit(1);
68  }
69 
71 
72  if (doSummary){
73  if (doDT){
74  dbe->setCurrentFolder(topFolder.str()+"/DT");
75  hLocalPositionDT=dbe->book2D("hLocalPositionDT","Local DT position (cm) absolute MEAN residuals;Sector;;cm", 14,1, 15,40,0,40);
76  hLocalAngleDT=dbe->book2D("hLocalAngleDT","Local DT angle (rad) absolute MEAN residuals;Sector;;rad", 14,1, 15,40,0,40);
77  hLocalPositionRmsDT=dbe->book2D("hLocalPositionRmsDT","Local DT position (cm) RMS residuals;Sector;;cm", 14,1, 15,40,0,40);
78  hLocalAngleRmsDT=dbe->book2D("hLocalAngleRmsDT","Local DT angle (rad) RMS residuals;Sector;;rad", 14,1, 15,40,0,40);
79 
80  hLocalXMeanDT=dbe->book1D("hLocalXMeanDT","Distribution of absolute MEAN Local X (cm) residuals for DT;<X> (cm);number of chambers",100,0,meanPositionRange);
81  hLocalXRmsDT=dbe->book1D("hLocalXRmsDT","Distribution of RMS Local X (cm) residuals for DT;X RMS (cm);number of chambers", 100,0,rmsPositionRange);
82  hLocalYMeanDT=dbe->book1D("hLocalYMeanDT","Distribution of absolute MEAN Local Y (cm) residuals for DT;<Y> (cm);number of chambers", 100,0,meanPositionRange);
83  hLocalYRmsDT=dbe->book1D("hLocalYRmsDT","Distribution of RMS Local Y (cm) residuals for DT;Y RMS (cm);number of chambers", 100,0,rmsPositionRange);
84 
85  hLocalPhiMeanDT=dbe->book1D("hLocalPhiMeanDT","Distribution of MEAN #phi (rad) residuals for DT;<#phi>(rad);number of chambers", 100,-meanAngleRange,meanAngleRange);
86  hLocalPhiRmsDT=dbe->book1D("hLocalPhiRmsDT","Distribution of RMS #phi (rad) residuals for DT;#phi RMS (rad);number of chambers", 100,0,rmsAngleRange);
87  hLocalThetaMeanDT=dbe->book1D("hLocalThetaMeanDT","Distribution of MEAN #theta (rad) residuals for DT;<#theta>(rad);number of chambers", 100,-meanAngleRange,meanAngleRange);
88  hLocalThetaRmsDT=dbe->book1D("hLocalThetaRmsDT","Distribution of RMS #theta (rad) residuals for DT;#theta RMS (rad);number of chambers",100,0,rmsAngleRange);
89  }
90 
91  if (doCSC){
92  dbe->setCurrentFolder(topFolder.str()+"/CSC");
93  hLocalPositionCSC=dbe->book2D("hLocalPositionCSC","Local CSC position (cm) absolute MEAN residuals;Sector;;cm",36,1,37,40,0,40);
94  hLocalAngleCSC=dbe->book2D("hLocalAngleCSC","Local CSC angle (rad) absolute MEAN residuals;Sector;;rad", 36,1,37,40,0,40);
95  hLocalPositionRmsCSC=dbe->book2D("hLocalPositionRmsCSC","Local CSC position (cm) RMS residuals;Sector;;cm", 36,1,37,40,0,40);
96  hLocalAngleRmsCSC=dbe->book2D("hLocalAngleRmsCSC","Local CSC angle (rad) RMS residuals;Sector;;rad", 36,1,37,40,0,40);
97 
98  hLocalXMeanCSC=dbe->book1D("hLocalXMeanCSC","Distribution of absolute MEAN Local X (cm) residuals for CSC;<X> (cm);number of chambers",100,0,meanPositionRange);
99  hLocalXRmsCSC=dbe->book1D("hLocalXRmsCSC","Distribution of RMS Local X (cm) residuals for CSC;X RMS (cm);number of chambers", 100,0,rmsPositionRange);
100  hLocalYMeanCSC=dbe->book1D("hLocalYMeanCSC","Distribution of absolute MEAN Local Y (cm) residuals for CSC;<Y> (cm);number of chambers", 100,0,meanPositionRange);
101  hLocalYRmsCSC=dbe->book1D("hLocalYRmsCSC","Distribution of RMS Local Y (cm) residuals for CSC;Y RMS (cm);number of chambers", 100,0,rmsPositionRange);
102 
103  hLocalPhiMeanCSC=dbe->book1D("hLocalPhiMeanCSC","Distribution of absolute MEAN #phi (rad) residuals for CSC;<#phi>(rad);number of chambers", 100,0,meanAngleRange);
104  hLocalPhiRmsCSC=dbe->book1D("hLocalPhiRmsCSC","Distribution of RMS #phi (rad) residuals for CSC;#phi RMS (rad);number of chambers", 100,0,rmsAngleRange);
105  hLocalThetaMeanCSC=dbe->book1D("hLocalThetaMeanCSC","Distribution of absolute MEAN #theta (rad) residuals for CSC;<#theta>(rad);number of chambers", 100,0,meanAngleRange);
106  hLocalThetaRmsCSC=dbe->book1D("hLocalThetaRmsCSC","Distribution of RMS #theta (rad) residuals for CSC;#theta RMS (rad);number of chambers",100,0,rmsAngleRange);
107  }
108  }
109 
110 
111  // Chamber individual histograms
112  // I need to create all of them even if they are empty to allow proper root merging
113 
114  // variables for histos ranges
115  double rangeX=0,rangeY=0;
116  std::string nameOfHistoLocalX,nameOfHistoLocalY,nameOfHistoLocalPhi,nameOfHistoLocalTheta;
117 
118  for (int station = -4; station<5; station++){
119 
120  //This piece of code calculates the range of the residuals
121  switch(abs(station)) {
122  case 1:
123  {rangeX = resLocalXRangeStation1; rangeY = resLocalYRangeStation1;}
124  break;
125  case 2:
126  {rangeX = resLocalXRangeStation2; rangeY = resLocalYRangeStation2;}
127  break;
128  case 3:
129  {rangeX = resLocalXRangeStation3; rangeY = resLocalYRangeStation3;}
130  break;
131  case 4:
132  {rangeX = resLocalXRangeStation4; rangeY = resLocalYRangeStation4;}
133  break;
134  default:
135  break;
136  }
137  if (doDT){
138  if(station>0){
139 
140  for(int wheel=-2;wheel<3;wheel++){
141 
142  for (int sector=1;sector<15;sector++){
143 
144  if(!((sector==13 || sector ==14) && station!=4)){
145 
146  std::stringstream Wheel; Wheel<<wheel;
147  std::stringstream Station; Station<<station;
148  std::stringstream Sector; Sector<<sector;
149 
150  nameOfHistoLocalX="ResidualLocalX_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
151  nameOfHistoLocalPhi= "ResidualLocalPhi_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
152  nameOfHistoLocalTheta= "ResidualLocalTheta_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
153  nameOfHistoLocalY= "ResidualLocalY_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
154 
156  "/DT/Wheel"+Wheel.str()+
157  "/Station"+Station.str()+
158  "/Sector"+Sector.str());
159 
160  //Create ME and push histos into their respective vectors
161 
162  MonitorElement *histoLocalX = dbe->book1D(nameOfHistoLocalX, nameOfHistoLocalX, nbins, -rangeX, rangeX);
163  unitsLocalX.push_back(histoLocalX);
164  MonitorElement *histoLocalPhi = dbe->book1D(nameOfHistoLocalPhi, nameOfHistoLocalPhi, nbins, -resPhiRange, resPhiRange);
165  unitsLocalPhi.push_back(histoLocalPhi);
166  MonitorElement *histoLocalTheta = dbe->book1D(nameOfHistoLocalTheta, nameOfHistoLocalTheta, nbins, -resThetaRange, resThetaRange);
167  unitsLocalTheta.push_back(histoLocalTheta);
168  MonitorElement *histoLocalY = dbe->book1D(nameOfHistoLocalY, nameOfHistoLocalY, nbins, -rangeY, rangeY);
169  unitsLocalY.push_back(histoLocalY);
170  }
171  }
172  }
173  } //station>0
174  }// doDT
175 
176  if (doCSC){
177  if(station!=0){
178 
179  for(int ring=1;ring<5;ring++){
180 
181  for(int chamber=1;chamber<37;chamber++){
182 
183  if( !( ((abs(station)==2 || abs(station)==3 || abs(station)==4) && ring==1 && chamber>18) ||
184  ((abs(station)==2 || abs(station)==3 || abs(station)==4) && ring>2)) ){
185  std::stringstream Ring; Ring<<ring;
186  std::stringstream Station; Station<<station;
187  std::stringstream Chamber; Chamber<<chamber;
188 
189  nameOfHistoLocalX="ResidualLocalX_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
190  nameOfHistoLocalPhi= "ResidualLocalPhi_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
191  nameOfHistoLocalTheta= "ResidualLocalTheta_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
192  nameOfHistoLocalY= "ResidualLocalY_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
193 
195  "/CSC/Station"+Station.str()+
196  "/Ring"+Ring.str()+
197  "/Chamber"+Chamber.str());
198 
199  //Create ME and push histos into their respective vectors
200 
201  MonitorElement *histoLocalX = dbe->book1D(nameOfHistoLocalX, nameOfHistoLocalX, nbins, -rangeX, rangeX);
202  unitsLocalX.push_back(histoLocalX);
203  MonitorElement *histoLocalPhi = dbe->book1D(nameOfHistoLocalPhi, nameOfHistoLocalPhi, nbins, -resPhiRange, resPhiRange);
204  unitsLocalPhi.push_back(histoLocalPhi);
205  MonitorElement *histoLocalTheta = dbe->book1D(nameOfHistoLocalTheta, nameOfHistoLocalTheta, nbins, -resThetaRange, resThetaRange);
206  unitsLocalTheta.push_back(histoLocalTheta);
207  MonitorElement *histoLocalY = dbe->book1D(nameOfHistoLocalY, nameOfHistoLocalY, nbins, -rangeY, rangeY);
208  unitsLocalY.push_back(histoLocalY);
209 
210  }
211  }
212  }
213  } // station!=0
214  }// doCSC
215 
216  } // loop on stations
217 }
218 
219 
221 
222  LogTrace(metname)<<"[MuonAlignment] Analysis of event # ";
223 
224  edm::ESHandle<MagneticField> theMGField;
225  iSetup.get<IdealMagneticFieldRecord>().get(theMGField);
226 
227  edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
228  iSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
229 
230  edm::ESHandle<Propagator> thePropagatorOpp;
231  iSetup.get<TrackingComponentsRecord>().get( "SmartPropagatorOpposite", thePropagatorOpp );
232 
233  edm::ESHandle<Propagator> thePropagatorAlo;
234  iSetup.get<TrackingComponentsRecord>().get( "SmartPropagator", thePropagatorAlo );
235 
236 // edm::ESHandle<Propagator> thePropagator;
237 // iSetup.get<TrackingComponentsRecord>().get( "SmartPropagatorAnyOpposite", thePropagator );
238 
239  // Get the RecoMuons collection from the event
241  iEvent.getByToken(theMuonCollectionLabel, muons);
242 
243  // Get the 4D DTSegments
245  iEvent.getByToken(theRecHits4DTagDT, all4DSegmentsDT);
247 
248  // Get the 4D CSCSegments
249  edm::Handle<CSCSegmentCollection> all4DSegmentsCSC;
250  iEvent.getByToken(theRecHits4DTagCSC, all4DSegmentsCSC);
252 
253  //Vectors used to perform the matching between Segments and hits from Track
254  intDVector indexCollectionDT;
255  intDVector indexCollectionCSC;
256 
257 
258 // thePropagator = new SteppingHelixPropagator(&*theMGField, alongMomentum);
259 
260  int countTracks = 0;
261  reco::TrackCollection::const_iterator muon;
262  for (muon = muons->begin(); muon != muons->end(); ++muon){
263 // if(muon->isGlobalMuon()){
264 // if(muon->isStandAloneMuon()){
265 
266  int countPoints = 0;
267 
268 // reco::TrackRef trackTR = muon->innerTrack();
269 // reco::TrackRef trackSA = muon->outerTrack();
270  //reco::TrackRef trackSA = muon;
271 
272  if(muon->recHitsSize()>(min1DTrackRecHitSize-1)) {
273 
274 // reco::TransientTrack tTrackTR( *trackTR, &*theMGField, theTrackingGeometry );
275  reco::TransientTrack tTrackSA(*muon,&*theMGField,theTrackingGeometry);
276 
277  // Adapted code for muonCosmics
278 
279  Double_t innerPerpSA = tTrackSA.innermostMeasurementState().globalPosition().perp();
280  Double_t outerPerpSA = tTrackSA.outermostMeasurementState().globalPosition().perp();
281 
283 // PropagationDirection propagationDir=alongMomentum;
284  const Propagator * thePropagator;
285 
286  // Define which kind of reco track is used
287  if ( (outerPerpSA-innerPerpSA) > 0 ) {
288 
289  trackRefitterType = "LHCLike";
290  innerTSOS = tTrackSA.innermostMeasurementState();
291  thePropagator = thePropagatorAlo.product();
292 // propagationDir = alongMomentum;
293 
294  }else {//if ((outerPerpSA-innerPerpSA) < 0 ) {
295 
296  trackRefitterType = "CosmicLike";
297  innerTSOS = tTrackSA.outermostMeasurementState();
298  thePropagator = thePropagatorOpp.product();
299 // propagationDir = oppositeToMomentum;
300 
301  }
302 
303  RecHitVector my4DTrack = this->doMatching(*muon, all4DSegmentsDT, all4DSegmentsCSC, &indexCollectionDT, &indexCollectionCSC, theTrackingGeometry);
304 
305 
306 //cut in number of segments
307  if(my4DTrack.size()>(min4DTrackSegmentSize-1) ){
308 
309 
310 // start propagation
311 // TrajectoryStateOnSurface innerTSOS = track.impactPointState();
312 // TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
313 
314  //If the state is valid
315  if(innerTSOS.isValid()) {
316 
317  //Loop over Associated segments
318  for(RecHitVector::iterator rechit = my4DTrack.begin(); rechit != my4DTrack.end(); ++rechit) {
319 
320  const GeomDet* geomDet = theTrackingGeometry->idToDet((*rechit)->geographicalId());
321 //Otherwise the propagator could throw an exception
322  const Plane* pDest = dynamic_cast<const Plane*>(&geomDet->surface());
323  const Cylinder* cDest = dynamic_cast<const Cylinder*>(&geomDet->surface());
324 
325  if(pDest != 0 || cDest != 0) {
326 
327 // Propagator *updatePropagator=thePropagator->clone();
328 // updatePropagator->setPropagationDirection(propagationDir);
329  TrajectoryStateOnSurface destiny = thePropagator->propagate(*(innerTSOS.freeState()), geomDet->surface());
330 
331  if(!destiny.isValid()|| !destiny.hasError()) continue;
332 
333  const long rawId= (*rechit)->geographicalId().rawId();
334  int position = -1;
335 
336  DetId myDet(rawId);
337  int det = myDet.subdetId();
338  int wheel=0,station=0,sector=0;
339  int endcap=0,ring=0,chamber=0;
340  bool goAhead = (det==1 && doDT) || (det==2 && doCSC);
341 
342  double residualLocalX=0,residualLocalPhi=0,residualLocalY=0,residualLocalTheta=0;
343 
344  // Fill generic histograms
345  //If it's a DT
346  if(det == 1 && doDT) {
347  DTChamberId myChamber(rawId);
348  wheel=myChamber.wheel();
349  station = myChamber.station();
350  sector=myChamber.sector();
351 
352  residualLocalX = (*rechit)->localPosition().x() -destiny.localPosition().x();
353 
354 
355  residualLocalPhi = atan2(((RecSegment *)(*rechit))->localDirection().z(),
356  ((RecSegment*)(*rechit))->localDirection().x()) - atan2(destiny.localDirection().z(), destiny.localDirection().x());
357  if(station!=4){
358 
359  residualLocalY = (*rechit)->localPosition().y() - destiny.localPosition().y();
360 
361 
362  residualLocalTheta = atan2(((RecSegment *)(*rechit))->localDirection().z(),
363  ((RecSegment*)(*rechit))->localDirection().y()) - atan2(destiny.localDirection().z(), destiny.localDirection().y());
364 
365  }
366 
367  }
368  else if (det==2 && doCSC){
369  CSCDetId myChamber(rawId);
370  endcap= myChamber.endcap();
371  station = myChamber.station();
372  if(endcap==2) station = -station;
373  ring = myChamber.ring();
374  chamber=myChamber.chamber();
375 
376  residualLocalX = (*rechit)->localPosition().x() -destiny.localPosition().x();
377 
378  residualLocalY = (*rechit)->localPosition().y() - destiny.localPosition().y();
379 
380  residualLocalPhi = atan2(((RecSegment *)(*rechit))->localDirection().y(),
381  ((RecSegment*)(*rechit))->localDirection().x()) - atan2(destiny.localDirection().y(), destiny.localDirection().x());
382 
383  residualLocalTheta = atan2(((RecSegment *)(*rechit))->localDirection().y(),
384  ((RecSegment*)(*rechit))->localDirection().z()) - atan2(destiny.localDirection().y(), destiny.localDirection().z());
385 
386  }
387  else{
388  residualLocalX=0,residualLocalPhi=0,residualLocalY=0,residualLocalTheta=0;
389  }
390 
391  // Fill individual chamber histograms
392 
393 
394  std::string nameOfHistoLocalX;
395 
396 
397  if(det==1 && doDT){ // DT
398  std::stringstream Wheel; Wheel<<wheel;
399  std::stringstream Station; Station<<station;
400  std::stringstream Sector; Sector<<sector;
401 
402  nameOfHistoLocalX="ResidualLocalX_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
403 
404  }
405  else if(det==2 && doCSC){ //CSC
406  std::stringstream Ring; Ring<<ring;
407  std::stringstream Station; Station<<station;
408  std::stringstream Chamber; Chamber<<chamber;
409 
410  nameOfHistoLocalX="ResidualLocalX_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
411  }
412 
413  if (goAhead){
414 
415  for(unsigned int i=0 ; i<unitsLocalX.size() ; i++)
416  {
417 
418  if( nameOfHistoLocalX==unitsLocalX[i]->getName()){
419  position=i; break;}
420  }
421 
422 
423  if(trackRefitterType == "CosmicLike") { //problem with angle convention in reverse extrapolation
424  residualLocalPhi += 3.1416;
425  residualLocalTheta +=3.1416;
426  }
427 
428  unitsLocalX.at(position)->Fill(residualLocalX);
429  unitsLocalPhi.at(position)->Fill(residualLocalPhi);
430 
431  if(det==1 && station!=4) {unitsLocalY.at(position)->Fill(residualLocalY);
432  unitsLocalTheta.at(position)->Fill(residualLocalTheta);}
433 
434  else if(det==2) {unitsLocalY.at(position)->Fill(residualLocalY);
435  unitsLocalTheta.at(position)->Fill(residualLocalTheta);}
436 
437 
438  countPoints++;
439  // if at least one point from this track is used, count this track
440  if (countPoints==1) countTracks++;
441  }
442 
443  innerTSOS = destiny;
444 
445  //delete thePropagator;
446 
447  }else {
448  edm::LogError("MuonAlignment") <<" Error!! Exception in propagator catched" << std::endl;
449  continue;
450  }
451 
452  } //loop over my4DTrack
453  } //TSOS was valid
454 
455  } // cut in at least 4 segments
456 
457  } //end cut in RecHitsSize>36
458  numberOfHits=numberOfHits+countPoints;
459 // } //Muon is GlobalMuon
460  } //loop over Muons
462 
463 // delete thePropagator;
464 
465 // else edm::LogError("MuonAlignment")<<"Error!! Specified MuonCollection "<< theMuonTrackCollectionLabel <<" is not present in event. ProductNotFound!!"<<std::endl;
466 }
467 
468 RecHitVector MuonAlignment::doMatching(const reco::Track &staTrack, edm::Handle<DTRecSegment4DCollection> &all4DSegmentsDT, edm::Handle<CSCSegmentCollection> &all4DSegmentsCSC, intDVector *indexCollectionDT, intDVector *indexCollectionCSC, edm::ESHandle<GlobalTrackingGeometry> &theTrackingGeometry) {
469 
472 
473  std::vector<int> positionDT;
474  std::vector<int> positionCSC;
475  RecHitVector my4DTrack;
476 
477  //Loop over the hits of the track
478  for(unsigned int counter = 0; counter != staTrack.recHitsSize()-1; counter++) {
479 
480  TrackingRecHitRef myRef = staTrack.recHit(counter);
481  const TrackingRecHit *rechit = myRef.get();
482  const GeomDet* geomDet = theTrackingGeometry->idToDet(rechit->geographicalId());
483 
484  //It's a DT Hit
485  if(geomDet->subDetector() == GeomDetEnumerators::DT) {
486 
487  //Take the layer associated to this hit
488  DTLayerId myLayer(rechit->geographicalId().rawId());
489 
490  int NumberOfDTSegment = 0;
491  //Loop over segments
492  for(segmentDT = all4DSegmentsDT->begin(); segmentDT != all4DSegmentsDT->end(); ++segmentDT) {
493 
494  //By default the chamber associated to this Segment is new
495  bool isNewChamber = true;
496 
497  //Loop over segments already included in the vector of segments in the actual track
498  for(std::vector<int>::iterator positionIt = positionDT.begin(); positionIt != positionDT.end(); positionIt++) {
499 
500  //If this segment has been used before isNewChamber = false
501  if(NumberOfDTSegment == *positionIt) isNewChamber = false;
502  }
503 
504  //Loop over vectors of segments associated to previous tracks
505  for(std::vector<std::vector<int> >::iterator collect = indexCollectionDT->begin(); collect != indexCollectionDT->end(); ++collect) {
506 
507  //Loop over segments associated to a track
508  for(std::vector<int>::iterator positionIt = (*collect).begin(); positionIt != (*collect).end(); positionIt++) {
509 
510  //If this segment was used in a previos track then isNewChamber = false
511  if(NumberOfDTSegment == *positionIt) isNewChamber = false;
512  }
513  }
514 
515  //If the chamber is new
516  if(isNewChamber) {
517 
518  DTChamberId myChamber((*segmentDT).geographicalId().rawId());
519  //If the layer of the hit belongs to the chamber of the 4D Segment
520  if(myLayer.wheel() == myChamber.wheel() && myLayer.station() == myChamber.station() && myLayer.sector() == myChamber.sector()) {
521 
522  //push position of the segment and tracking rechit
523  positionDT.push_back(NumberOfDTSegment);
524  my4DTrack.push_back((TrackingRecHit *) &(*segmentDT));
525  }
526  }
527  NumberOfDTSegment++;
528  }
529  //In case is a CSC
530  } else if (geomDet->subDetector() == GeomDetEnumerators::CSC) {
531 
532  //Take the layer associated to this hit
533  CSCDetId myLayer(rechit->geographicalId().rawId());
534 
535  int NumberOfCSCSegment = 0;
536  //Loop over 4Dsegments
537  for(segmentCSC = all4DSegmentsCSC->begin(); segmentCSC != all4DSegmentsCSC->end(); segmentCSC++) {
538 
539  //By default the chamber associated to the segment is new
540  bool isNewChamber = true;
541 
542  //Loop over segments in the current track
543  for(std::vector<int>::iterator positionIt = positionCSC.begin(); positionIt != positionCSC.end(); positionIt++) {
544 
545  //If this segment has been used then newchamber = false
546  if(NumberOfCSCSegment == *positionIt) isNewChamber = false;
547  }
548  //Loop over vectors of segments in previous tracks
549  for(std::vector<std::vector<int> >::iterator collect = indexCollectionCSC->begin(); collect != indexCollectionCSC->end(); ++collect) {
550 
551  //Loop over segments in a track
552  for(std::vector<int>::iterator positionIt = (*collect).begin(); positionIt != (*collect).end(); positionIt++) {
553 
554  //If the segment was used in a previous track isNewChamber = false
555  if(NumberOfCSCSegment == *positionIt) isNewChamber = false;
556  }
557  }
558  //If the chamber is new
559  if(isNewChamber) {
560 
561  CSCDetId myChamber((*segmentCSC).geographicalId().rawId());
562  //If the chambers are the same
563  if(myLayer.chamberId() == myChamber.chamberId()) {
564  //push
565  positionCSC.push_back(NumberOfCSCSegment);
566  my4DTrack.push_back((TrackingRecHit *) &(*segmentCSC));
567  }
568  }
569  NumberOfCSCSegment++;
570  }
571  }
572  }
573 
574  indexCollectionDT->push_back(positionDT);
575  indexCollectionCSC->push_back(positionCSC);
576 
577  if ( trackRefitterType == "CosmicLike") {
578 
579  std::reverse(my4DTrack.begin(),my4DTrack.end());
580 
581  }
582  return my4DTrack;
583 
584 
585 }
586 
587 
588 
590 
591 
592  LogTrace(metname)<<"[MuonAlignment] Saving the histos";
593  bool outputMEsInRootFile = parameters.getParameter<bool>("OutputMEsInRootFile");
595 
596  edm::LogInfo("MuonAlignment") << "Number of Tracks considered for residuals: " << numberOfTracks << std::endl << std::endl;
597  edm::LogInfo("MuonAlignment") << "Number of Hits considered for residuals: " << numberOfHits << std::endl << std::endl;
598 
599  if (doSummary){
600  char binLabel[15];
601 
602  // check if ME still there (and not killed by MEtoEDM for memory saving)
603  if( dbe )
604  {
605  // check existence of first histo in the list
606  if (! dbe->get(topFolder.str()+"/DT/hLocalPositionDT")) return;
607  }
608  else
609  return;
610 
611 
612  for(unsigned int i=0 ; i<unitsLocalX.size() ; i++)
613  {
614 
615  if(unitsLocalX[i]->getEntries()!=0){
616 
617  TString nameHistoLocalX=unitsLocalX[i]->getName();
618 
619  TString nameHistoLocalPhi=unitsLocalPhi[i]->getName();
620 
621  TString nameHistoLocalTheta=unitsLocalTheta[i]->getName();
622 
623  TString nameHistoLocalY=unitsLocalY[i]->getName();
624 
625 
626  if (nameHistoLocalX.Contains("MB") ) // HistoLocalX DT
627  {
628  int wheel, station, sector;
629 
630  sscanf(nameHistoLocalX, "ResidualLocalX_W%dMB%1dS%d",&wheel,&station,&sector);
631 
632  Int_t nstation=station - 1;
633  Int_t nwheel=wheel+2;
634 
635  Double_t Mean=unitsLocalX[i]->getMean();
636  Double_t Error=unitsLocalX[i]->getMeanError();
637 
638  Int_t ybin=1+nwheel*8+nstation*2;
639  hLocalPositionDT->setBinContent(sector,ybin,fabs(Mean));
640  snprintf(binLabel, sizeof(binLabel), "MB%d/%d_X",wheel, station );
641  hLocalPositionDT->setBinLabel(ybin,binLabel,2);
642  hLocalPositionRmsDT->setBinContent(sector,ybin,Error);
643  hLocalPositionRmsDT->setBinLabel(ybin,binLabel,2);
644 
645  hLocalXMeanDT->Fill(fabs(Mean));
646  hLocalXRmsDT->Fill(Error);
647  }
648 
649  if (nameHistoLocalX.Contains("ME")) // HistoLocalX CSC
650  {
651  int station, ring, chamber;
652 
653  sscanf(nameHistoLocalX, "ResidualLocalX_ME%dR%1dC%d",&station,&ring,&chamber);
654 
655  Double_t Mean=unitsLocalX[i]->getMean();
656  Double_t Error=unitsLocalX[i]->getMeanError();
657 
658  Int_t ybin=abs(station)*2+ring;
659  if(abs(station)==1) ybin=ring;
660  if (station>0) ybin=ybin+10;
661  else ybin = 11 -ybin;
662  ybin=2*ybin-1;
663  hLocalPositionCSC->setBinContent(chamber,ybin,fabs(Mean));
664  snprintf(binLabel, sizeof(binLabel), "ME%d/%d_X", station,ring );
665  hLocalPositionCSC->setBinLabel(ybin,binLabel,2);
666  hLocalPositionRmsCSC->setBinContent(chamber,ybin,Error);
667  hLocalPositionRmsCSC->setBinLabel(ybin,binLabel,2);
668 
669  hLocalXMeanCSC->Fill(fabs(Mean));
670  hLocalXRmsCSC->Fill(Error);
671  }
672 
673  if (nameHistoLocalTheta.Contains("MB")) // HistoLocalTheta DT
674  {
675 
676  int wheel, station, sector;
677 
678  sscanf(nameHistoLocalTheta, "ResidualLocalTheta_W%dMB%1dS%d",&wheel,&station,&sector);
679 
680  if(station != 4){
681  Int_t nstation=station - 1;
682  Int_t nwheel=wheel+2;
683 
684  Double_t Mean=unitsLocalTheta[i]->getMean();
685  Double_t Error=unitsLocalTheta[i]->getMeanError();
686 
687  Int_t ybin=2+nwheel*8+nstation*2;
688  hLocalAngleDT->setBinContent(sector,ybin,fabs(Mean));
689  snprintf(binLabel, sizeof(binLabel), "MB%d/%d_#theta",wheel,station );
690  hLocalAngleDT->setBinLabel(ybin,binLabel,2);
691  hLocalAngleRmsDT->setBinContent(sector,ybin,Error);
692  hLocalAngleRmsDT->setBinLabel(ybin,binLabel,2);
693 
694  hLocalThetaMeanDT->Fill(fabs(Mean));
695  hLocalThetaRmsDT->Fill(Error);
696  }
697  }
698 
699  if (nameHistoLocalPhi.Contains("MB")) // HistoLocalPhi DT
700  {
701 
702  int wheel, station, sector;
703 
704  sscanf(nameHistoLocalPhi, "ResidualLocalPhi_W%dMB%1dS%d",&wheel,&station,&sector);
705 
706  Int_t nstation=station - 1;
707  Int_t nwheel=wheel+2;
708 
709  Double_t Mean=unitsLocalPhi[i]->getMean();
710  Double_t Error=unitsLocalPhi[i]->getMeanError();
711 
712  Int_t ybin=1+nwheel*8+nstation*2;
713  hLocalAngleDT->setBinContent(sector,ybin,fabs(Mean));
714  snprintf(binLabel, sizeof(binLabel), "MB%d/%d_#phi", wheel,station );
715  hLocalAngleDT->setBinLabel(ybin,binLabel,2);
716  hLocalAngleRmsDT->setBinContent(sector,ybin,Error);
717  hLocalAngleRmsDT->setBinLabel(ybin,binLabel,2);
718 
719  hLocalPhiMeanDT->Fill(fabs(Mean));
720  hLocalPhiRmsDT->Fill(Error);
721  }
722 
723  if (nameHistoLocalPhi.Contains("ME")) // HistoLocalPhi CSC
724  {
725 
726  int station, ring, chamber;
727 
728  sscanf(nameHistoLocalPhi, "ResidualLocalPhi_ME%dR%1dC%d",&station,&ring,&chamber);
729 
730  Double_t Mean=unitsLocalPhi[i]->getMean();
731  Double_t Error=unitsLocalPhi[i]->getMeanError();
732 
733  Int_t ybin=abs(station)*2+ring;
734  if(abs(station)==1) ybin=ring;
735  if (station>0) ybin=ybin+10;
736  else ybin = 11 -ybin;
737  ybin=2*ybin-1;
738  hLocalAngleCSC->setBinContent(chamber,ybin,fabs(Mean));
739  snprintf(binLabel, sizeof(binLabel), "ME%d/%d_#phi", station,ring );
740  hLocalAngleCSC->setBinLabel(ybin,binLabel,2);
741  hLocalAngleRmsCSC->setBinContent(chamber,ybin,Error);
742  hLocalAngleRmsCSC->setBinLabel(ybin,binLabel,2);
743 
744  hLocalPhiMeanCSC->Fill(fabs(Mean));
745  hLocalPhiRmsCSC->Fill(Error);
746  }
747 
748  if (nameHistoLocalTheta.Contains("ME")) // HistoLocalTheta CSC
749  {
750 
751  int station, ring, chamber;
752 
753  sscanf(nameHistoLocalTheta, "ResidualLocalTheta_ME%dR%1dC%d",&station,&ring,&chamber);
754 
755  Double_t Mean=unitsLocalTheta[i]->getMean();
756  Double_t Error=unitsLocalTheta[i]->getMeanError();
757 
758  Int_t ybin=abs(station)*2+ring;
759  if(abs(station)==1) ybin=ring;
760  if (station>0) ybin=ybin+10;
761  else ybin = 11 -ybin;
762  ybin=2*ybin;
763  hLocalAngleCSC->setBinContent(chamber,ybin,fabs(Mean));
764  snprintf(binLabel, sizeof(binLabel), "ME%d/%d_#theta", station,ring );
765  hLocalAngleCSC->setBinLabel(ybin,binLabel,2);
766  hLocalAngleRmsCSC->setBinContent(chamber,ybin,Error);
767  hLocalAngleRmsCSC->setBinLabel(ybin,binLabel,2);
768 
769  hLocalThetaMeanCSC->Fill(fabs(Mean));
770  hLocalThetaRmsCSC->Fill(Error);
771 
772  }
773 
774  if (nameHistoLocalY.Contains("MB")) // HistoLocalY DT
775  {
776 
777  int wheel, station, sector;
778 
779  sscanf(nameHistoLocalY, "ResidualLocalY_W%dMB%1dS%d",&wheel,&station,&sector);
780 
781  if(station!=4){
782  Int_t nstation=station - 1;
783  Int_t nwheel=wheel+2;
784 
785  Double_t Mean=unitsLocalY[i]->getMean();
786  Double_t Error=unitsLocalY[i]->getMeanError();
787 
788  Int_t ybin=2+nwheel*8+nstation*2;
789  hLocalPositionDT->setBinContent(sector,ybin,fabs(Mean));
790  snprintf(binLabel, sizeof(binLabel), "MB%d/%d_Y", wheel,station );
791  hLocalPositionDT->setBinLabel(ybin,binLabel,2);
792  hLocalPositionRmsDT->setBinContent(sector,ybin,Error);
793  hLocalPositionRmsDT->setBinLabel(ybin,binLabel,2);
794 
795  hLocalYMeanDT->Fill(fabs(Mean));
796  hLocalYRmsDT->Fill(Error);
797  }
798  }
799 
800  if (nameHistoLocalY.Contains("ME")) // HistoLocalY CSC
801  {
802 
803  int station, ring, chamber;
804 
805  sscanf(nameHistoLocalY, "ResidualLocalY_ME%dR%1dC%d",&station,&ring,&chamber);
806 
807  Double_t Mean=unitsLocalY[i]->getMean();
808  Double_t Error=unitsLocalY[i]->getMeanError();
809 
810  Int_t ybin=abs(station)*2+ring;
811  if(abs(station)==1) ybin=ring;
812  if (station>0) ybin=ybin+10;
813  else ybin = 11 -ybin;
814  ybin=2*ybin;
815  hLocalPositionCSC->setBinContent(chamber,ybin,fabs(Mean));
816  snprintf(binLabel, sizeof(binLabel), "ME%d/%d_Y", station,ring );
817  hLocalPositionCSC->setBinLabel(ybin,binLabel,2);
818  hLocalPositionRmsCSC->setBinContent(chamber,ybin,Error);
819  hLocalPositionRmsCSC->setBinLabel(ybin,binLabel,2);
820 
821  hLocalYMeanCSC->Fill(fabs(Mean));
822  hLocalYRmsCSC->Fill(Error);
823  }
824  } // check in # entries
825  } // loop on vector of histos
826  } //doSummary
827 
828  if(outputMEsInRootFile){
829 // dbe->showDirStructure();
830  dbe->save(outputFileName);
831  }
832 
833 
834 }
835 
MonitorElement * hLocalPhiRmsDT
Definition: MuonAlignment.h:99
edm::EDGetTokenT< CSCSegmentCollection > theRecHits4DTagCSC
int chamber() const
Definition: CSCDetId.h:81
std::vector< MonitorElement * > unitsLocalX
MonitorElement * hLocalPositionRmsDT
Definition: MuonAlignment.h:90
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
virtual FreeTrajectoryState propagate(const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const final
Definition: Propagator.h:119
int i
Definition: DBlmapReader.cc:9
MonitorElement * hLocalThetaMeanCSC
MonitorElement * hLocalYMeanDT
Definition: MuonAlignment.h:96
void setBinContent(int binx, double content)
set content of bin (1-D)
RecHitVector doMatching(const reco::Track &, edm::Handle< DTRecSegment4DCollection > &, edm::Handle< CSCSegmentCollection > &, intDVector *, intDVector *, edm::ESHandle< GlobalTrackingGeometry > &)
MonitorElement * hLocalPhiRmsCSC
double resLocalYRangeStation2
T perp() const
Definition: PV3DBase.h:72
DQMStore * dbe
Definition: MuonAlignment.h:87
MonitorElement * hLocalYRmsCSC
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:994
MonitorElement * hLocalPositionCSC
double resLocalXRangeStation1
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
MonitorElement * hLocalXMeanDT
Definition: MuonAlignment.h:94
size_t recHitsSize() const
Get number of RecHits. (Warning, this includes invalid hits, which are not physical hits)...
Definition: Track.h:119
LocalVector localDirection() const
double resLocalYRangeStation1
std::string trackRefitterType
double resPhiRange
MonitorElement * hLocalPhiMeanDT
Definition: MuonAlignment.h:98
T y() const
Definition: PV3DBase.h:63
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
GlobalPoint globalPosition() const
void analyze(const edm::Event &, const edm::EventSetup &)
Get the analysis.
edm::EDGetTokenT< reco::TrackCollection > theMuonCollectionLabel
void beginJob()
Inizialize parameters for histo binning.
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
double meanAngleRange
Definition: Plane.h:17
std::vector< TrackingRecHit * > RecHitVector
MuonAlignment(const edm::EventSetup &iSetup)
double resThetaRange
TrajectoryStateOnSurface innermostMeasurementState() const
MonitorElement * hLocalAngleRmsCSC
void Fill(long long x)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
int endcap() const
Definition: CSCDetId.h:106
edm::ParameterSet parameters
int iEvent
Definition: GenABIO.cc:230
std::string MEFolderName
MonitorElement * hLocalAngleRmsDT
Definition: MuonAlignment.h:92
double resLocalXRangeStation2
FreeTrajectoryState const * freeState(bool withErrors=true) const
T z() const
Definition: PV3DBase.h:64
std::string metname
MonitorElement * hLocalAngleDT
Definition: MuonAlignment.h:91
MonitorElement * hLocalThetaRmsDT
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< MonitorElement * > unitsLocalTheta
MonitorElement * hLocalYRmsDT
Definition: MuonAlignment.h:97
MonitorElement * hLocalThetaRmsCSC
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:244
TrajectoryStateOnSurface outermostMeasurementState() const
MonitorElement * hLocalXMeanCSC
MonitorElement * get(const std::string &path) const
get ME from full pathname (e.g. &quot;my/long/dir/my_histo&quot;)
Definition: DQMStore.cc:1748
MonitorElement * hLocalYMeanCSC
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
#define LogTrace(id)
unsigned int nbins
std::vector< MonitorElement * > unitsLocalPhi
double rmsPositionRange
int ring() const
Definition: CSCDetId.h:88
Definition: DetId.h:18
std::vector< std::vector< int > > intDVector
double rmsAngleRange
std::vector< MonitorElement * > unitsLocalY
MonitorElement * hLocalPhiMeanCSC
edm::EDGetTokenT< DTRecSegment4DCollection > theRecHits4DTagDT
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
MonitorElement * hLocalAngleCSC
double resLocalYRangeStation3
MonitorElement * hLocalXRmsDT
Definition: MuonAlignment.h:95
unsigned int min1DTrackRecHitSize
void endJob(void)
Save the histos.
double resLocalYRangeStation4
tuple muons
Definition: patZpeak.py:38
double meanPositionRange
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, const uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE", const bool resetMEsAfterWriting=false)
Definition: DQMStore.cc:2602
int sector() const
Definition: DTChamberId.h:61
MonitorElement * hLocalPositionRmsCSC
static std::atomic< unsigned int > counter
static int position[264][3]
Definition: ReadPGInfo.cc:509
TrackingRecHitRef recHit(size_t i) const
Get i-th hit on the track.
Definition: Track.h:114
int station() const
Definition: CSCDetId.h:99
DetId geographicalId() const
MonitorElement * hLocalXRmsCSC
int station() const
Return the station number.
Definition: DTChamberId.h:51
unsigned int min4DTrackSegmentSize
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:1122
virtual SubDetector subDetector() const
Which subdetector.
Definition: GeomDet.cc:49
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:45
std::stringstream topFolder
T x() const
Definition: PV3DBase.h:62
MonitorElement * hLocalThetaMeanDT
MonitorElement * hLocalPositionDT
Definition: MuonAlignment.h:89
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:707
double resLocalXRangeStation4
double resLocalXRangeStation3