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