CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
MuonAlignmentFromReference.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: MuonAlignmentAlgorithms
4 // Class: MuonAlignmentFromReference
5 //
13 //
14 // Original Author: Jim Pivarski,,,
15 // Created: Sat Jan 24 16:20:28 CST 2009
16 // $Id: MuonAlignmentFromReference.cc,v 1.39 2011/10/13 00:03:12 khotilov Exp $
17 
18 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmBase.h"
19 
20 #include "FWCore/Framework/interface/Event.h"
21 #include "FWCore/Framework/interface/ESHandle.h"
22 #include "FWCore/MessageLogger/interface/MessageLogger.h"
23 #include "FWCore/ParameterSet/interface/ParameterSet.h"
24 #include "FWCore/Utilities/interface/InputTag.h"
25 #include "FWCore/ServiceRegistry/interface/Service.h"
26 #include "CommonTools/UtilAlgos/interface/TFileService.h"
27 
28 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
29 #include "Geometry/Records/interface/MuonGeometryRecord.h"
30 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
31 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
32 
33 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
34 #include "MagneticField/Engine/interface/MagneticField.h"
35 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
36 #include "TrackingTools/Records/interface/DetIdAssociatorRecord.h"
37 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
38 
39 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
40 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
41 #include "DataFormats/MuonDetId/interface/DTSuperLayerId.h"
42 #include "DataFormats/TrackReco/interface/Track.h"
43 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
44 
45 #include "TrackingTools/PatternTools/interface/Trajectory.h"
46 
47 #include "Alignment/CommonAlignment/interface/Alignable.h"
48 #include "Alignment/CommonAlignment/interface/AlignableDetUnit.h"
49 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
50 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
51 #include "Alignment/MuonAlignment/interface/AlignableDTSuperLayer.h"
52 #include "Alignment/MuonAlignment/interface/AlignableDTChamber.h"
53 #include "Alignment/MuonAlignment/interface/AlignableDTStation.h"
54 #include "Alignment/MuonAlignment/interface/AlignableDTWheel.h"
55 #include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
56 #include "Alignment/MuonAlignment/interface/AlignableCSCRing.h"
57 #include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"
58 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
59 
60 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
61 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
62 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals5DOFFitter.h"
63 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFrphiFitter.h"
64 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h"
65 
66 #include "TFile.h"
67 #include "TTree.h"
68 #include "TStopwatch.h"
69 
70 #include <map>
71 #include <sstream>
72 #include <fstream>
73 
74 class MuonAlignmentFromReference : public AlignmentAlgorithmBase {
75 public:
76  MuonAlignmentFromReference(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC);
77  ~MuonAlignmentFromReference() override;
78 
79  void initialize(const edm::EventSetup& iSetup,
80  AlignableTracker* alignableTracker,
81  AlignableMuon* alignableMuon,
82  AlignableExtras* extras,
83  AlignmentParameterStore* alignmentParameterStore) override;
84 
85  void startNewLoop() override{};
86 
87  void run(const edm::EventSetup& iSetup, const EventInfo& eventInfo) override;
88 
89  void processMuonResidualsFromTrack(MuonResidualsFromTrack& mrft);
90 
91  void terminate(const edm::EventSetup& iSetup) override;
92 
93 private:
94  bool numeric(std::string s);
95  int number(std::string s);
96  std::string chamberPrettyNameFromId(unsigned int idx);
97 
98  void parseReference(align::Alignables& reference,
99  const align::Alignables& all_DT_chambers,
100  const align::Alignables& all_CSC_chambers);
101 
102  void fitAndAlign();
103  void readTmpFiles();
104  void writeTmpFiles();
105 
106  void selectResidualsPeaks();
107  void correctBField();
108  void fiducialCuts();
109  void eraseNotSelectedResiduals();
110 
111  void fillNtuple();
112 
113  // tokens
114  const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> m_cscGeometryToken;
115  const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_globTackingToken;
116  const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagFieldToken;
117  const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_propToken;
118  const edm::ESGetToken<DetIdAssociator, DetIdAssociatorRecord> m_DetIdToken;
119  const MuonResidualsFromTrack::BuilderToken m_builderToken;
120 
121  // configutarion paramenters:
122  edm::InputTag m_muonCollectionTag;
123  std::vector<std::string> m_reference;
124  double m_minTrackPt;
125  double m_maxTrackPt;
126  double m_minTrackP;
127  double m_maxTrackP;
128  double m_maxDxy;
129  int m_minTrackerHits;
130  double m_maxTrackerRedChi2;
131  bool m_allowTIDTEC;
132  int m_minNCrossedChambers;
133  int m_minDT13Hits;
134  int m_minDT2Hits;
135  int m_minCSCHits;
136  std::string m_writeTemporaryFile;
137  std::vector<std::string> m_readTemporaryFiles;
138  bool m_doAlignment;
139  int m_strategy;
140  std::string m_residualsModel;
141  int m_minAlignmentHits;
142  bool m_twoBin;
143  bool m_combineME11;
144  bool m_weightAlignment;
145  std::string m_reportFileName;
146  double m_maxResSlopeY;
147  bool m_createNtuple;
148  double m_peakNSigma;
149  int m_BFieldCorrection;
150  bool m_doDT;
151  bool m_doCSC;
152  std::string m_useResiduals;
153 
154  // utility objects
155  AlignableNavigator* m_alignableNavigator;
156  AlignmentParameterStore* m_alignmentParameterStore;
157  align::Alignables m_alignables;
158  std::map<Alignable*, Alignable*> m_me11map;
159  std::map<Alignable*, MuonResidualsTwoBin*> m_fitters;
160  std::vector<unsigned int> m_indexes;
161  std::map<unsigned int, MuonResidualsTwoBin*> m_fitterOrder;
162 
163  // counters
164  long m_counter_events;
165  long m_counter_tracks;
166  long m_counter_trackmomentum;
167  long m_counter_trackdxy;
168  long m_counter_trackerhits;
169  long m_counter_trackerchi2;
170  long m_counter_trackertidtec;
171  long m_counter_minchambers;
172  long m_counter_totchambers;
173  long m_counter_station123;
174  long m_counter_station123valid;
175  long m_counter_station123dt13hits;
176  long m_counter_station123dt2hits;
177  long m_counter_station123aligning;
178  long m_counter_station4;
179  long m_counter_station4valid;
180  long m_counter_station4hits;
181  long m_counter_station4aligning;
182  long m_counter_csc;
183  long m_counter_cscvalid;
184  long m_counter_cschits;
185  long m_counter_cscaligning;
186  long m_counter_resslopey;
187 
188  // debug ntuple
189  void bookNtuple();
190  TTree* m_ttree;
191  MuonResidualsFitter::MuonAlignmentTreeRow m_tree_row;
192 
193  bool m_debug;
194 };
195 
196 MuonAlignmentFromReference::MuonAlignmentFromReference(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
197  : AlignmentAlgorithmBase(cfg, iC),
198  m_cscGeometryToken(iC.esConsumes<edm::Transition::BeginRun>()),
199  m_globTackingToken(iC.esConsumes()),
200  m_MagFieldToken(iC.esConsumes()),
201  m_propToken(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
202  m_DetIdToken(iC.esConsumes(edm::ESInputTag("", "MuonDetIdAssociator"))),
203  m_builderToken(iC.esConsumes(MuonResidualsFromTrack::builderESInputTag())),
204  m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag")),
205  m_reference(cfg.getParameter<std::vector<std::string> >("reference")),
206  m_minTrackPt(cfg.getParameter<double>("minTrackPt")),
207  m_maxTrackPt(cfg.getParameter<double>("maxTrackPt")),
208  m_minTrackP(cfg.getParameter<double>("minTrackP")),
209  m_maxTrackP(cfg.getParameter<double>("maxTrackP")),
210  m_maxDxy(cfg.getParameter<double>("maxDxy")),
211  m_minTrackerHits(cfg.getParameter<int>("minTrackerHits")),
212  m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2")),
213  m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC")),
214  m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers")),
215  m_minDT13Hits(cfg.getParameter<int>("minDT13Hits")),
216  m_minDT2Hits(cfg.getParameter<int>("minDT2Hits")),
217  m_minCSCHits(cfg.getParameter<int>("minCSCHits")),
218  m_writeTemporaryFile(cfg.getParameter<std::string>("writeTemporaryFile")),
219  m_readTemporaryFiles(cfg.getParameter<std::vector<std::string> >("readTemporaryFiles")),
220  m_doAlignment(cfg.getParameter<bool>("doAlignment")),
221  m_strategy(cfg.getParameter<int>("strategy")),
222  m_residualsModel(cfg.getParameter<std::string>("residualsModel")),
223  m_minAlignmentHits(cfg.getParameter<int>("minAlignmentHits")),
224  m_twoBin(cfg.getParameter<bool>("twoBin")),
225  m_combineME11(cfg.getParameter<bool>("combineME11")),
226  m_weightAlignment(cfg.getParameter<bool>("weightAlignment")),
227  m_reportFileName(cfg.getParameter<std::string>("reportFileName")),
228  m_maxResSlopeY(cfg.getParameter<double>("maxResSlopeY")),
229  m_createNtuple(cfg.getParameter<bool>("createNtuple")),
230  m_peakNSigma(cfg.getParameter<double>("peakNSigma")),
231  m_BFieldCorrection(cfg.getParameter<int>("bFieldCorrection")),
232  m_doDT(cfg.getParameter<bool>("doDT")),
233  m_doCSC(cfg.getParameter<bool>("doCSC")),
234  m_useResiduals(cfg.getParameter<std::string>("useResiduals")) {
235  // alignment requires a TFile to provide plots to check the fit output
236  // just filling the residuals lists does not
237  // but we don't want to wait until the end of the job to find out that the TFile is missing
238  if (m_doAlignment || m_createNtuple) {
239  edm::Service<TFileService> fs;
240  TFile& tfile = fs->file();
241  tfile.ls();
242  }
243 
244  m_ttree = nullptr;
245  if (m_createNtuple)
246  bookNtuple();
247 
248  m_counter_events = 0;
249  m_counter_tracks = 0;
250  m_counter_trackmomentum = 0;
251  m_counter_trackdxy = 0;
252  m_counter_trackerhits = 0;
253  m_counter_trackerchi2 = 0;
254  m_counter_trackertidtec = 0;
255  m_counter_minchambers = 0;
256  m_counter_totchambers = 0;
257  m_counter_station123 = 0;
258  m_counter_station123valid = 0;
259  m_counter_station123dt13hits = 0;
260  m_counter_station123dt2hits = 0;
261  m_counter_station123aligning = 0;
262  m_counter_station4 = 0;
263  m_counter_station4valid = 0;
264  m_counter_station4hits = 0;
265  m_counter_station4aligning = 0;
266  m_counter_csc = 0;
267  m_counter_cscvalid = 0;
268  m_counter_cschits = 0;
269  m_counter_cscaligning = 0;
270  m_counter_resslopey = 0;
271 
272  m_debug = false;
273 }
274 
275 MuonAlignmentFromReference::~MuonAlignmentFromReference() { delete m_alignableNavigator; }
276 
277 void MuonAlignmentFromReference::bookNtuple() {
278  edm::Service<TFileService> fs;
279  m_ttree = fs->make<TTree>("mual_ttree", "mual_ttree");
280  m_ttree->Branch("is_plus", &m_tree_row.is_plus, "is_plus/O");
281  m_ttree->Branch("is_dt", &m_tree_row.is_dt, "is_dt/O");
282  m_ttree->Branch("station", &m_tree_row.station, "station/b");
283  m_ttree->Branch("ring_wheel", &m_tree_row.ring_wheel, "ring_wheel/B");
284  m_ttree->Branch("sector", &m_tree_row.sector, "sector/b");
285  m_ttree->Branch("res_x", &m_tree_row.res_x, "res_x/F");
286  m_ttree->Branch("res_y", &m_tree_row.res_y, "res_y/F");
287  m_ttree->Branch("res_slope_x", &m_tree_row.res_slope_x, "res_slope_x/F");
288  m_ttree->Branch("res_slope_y", &m_tree_row.res_slope_y, "res_slope_y/F");
289  m_ttree->Branch("pos_x", &m_tree_row.pos_x, "pos_x/F");
290  m_ttree->Branch("pos_y", &m_tree_row.pos_y, "pos_y/F");
291  m_ttree->Branch("angle_x", &m_tree_row.angle_x, "angle_x/F");
292  m_ttree->Branch("angle_y", &m_tree_row.angle_y, "angle_y/F");
293  m_ttree->Branch("pz", &m_tree_row.pz, "pz/F");
294  m_ttree->Branch("pt", &m_tree_row.pt, "pt/F");
295  m_ttree->Branch("q", &m_tree_row.q, "q/B");
296  m_ttree->Branch("select", &m_tree_row.select, "select/O");
297  //m_ttree->Branch("",&m_tree_row.,"/");
298 }
299 
300 bool MuonAlignmentFromReference::numeric(std::string s) { return s.length() == 1 && std::isdigit(s[0]); }
301 
302 int MuonAlignmentFromReference::number(std::string s) {
303  if (!numeric(s))
304  assert(false);
305  return atoi(s.c_str());
306 }
307 
308 void MuonAlignmentFromReference::initialize(const edm::EventSetup& iSetup,
309  AlignableTracker* alignableTracker,
310  AlignableMuon* alignableMuon,
311  AlignableExtras* extras,
312  AlignmentParameterStore* alignmentParameterStore) {
313  if (alignableMuon == nullptr)
314  throw cms::Exception("MuonAlignmentFromReference") << "doMuon must be set to True" << std::endl;
315 
316  m_alignableNavigator = new AlignableNavigator(alignableMuon);
317  m_alignmentParameterStore = alignmentParameterStore;
318  m_alignables = m_alignmentParameterStore->alignables();
319 
320  int residualsModel;
321  if (m_residualsModel == std::string("pureGaussian"))
322  residualsModel = MuonResidualsFitter::kPureGaussian;
323  else if (m_residualsModel == std::string("pureGaussian2D"))
324  residualsModel = MuonResidualsFitter::kPureGaussian2D;
325  else if (m_residualsModel == std::string("powerLawTails"))
326  residualsModel = MuonResidualsFitter::kPowerLawTails;
327  else if (m_residualsModel == std::string("ROOTVoigt"))
328  residualsModel = MuonResidualsFitter::kROOTVoigt;
329  else if (m_residualsModel == std::string("GaussPowerTails"))
330  residualsModel = MuonResidualsFitter::kGaussPowerTails;
331  else
332  throw cms::Exception("MuonAlignmentFromReference")
333  << "unrecognized residualsModel: \"" << m_residualsModel << "\"" << std::endl;
334 
335  int useResiduals;
336  if (m_useResiduals == std::string("1111"))
337  useResiduals = MuonResidualsFitter::k1111;
338  else if (m_useResiduals == std::string("1110"))
339  useResiduals = MuonResidualsFitter::k1110;
340  else if (m_useResiduals == std::string("1100"))
341  useResiduals = MuonResidualsFitter::k1100;
342  else if (m_useResiduals == std::string("1000"))
343  useResiduals = MuonResidualsFitter::k1000;
344  else if (m_useResiduals == std::string("1010"))
345  useResiduals = MuonResidualsFitter::k1010;
346  else if (m_useResiduals == std::string("0010"))
347  useResiduals = MuonResidualsFitter::k0010;
348  else
349  throw cms::Exception("MuonAlignmentFromReference")
350  << "unrecognized useResiduals: \"" << m_useResiduals << "\"" << std::endl;
351 
352  const CSCGeometry* cscGeometry = &iSetup.getData(m_cscGeometryToken);
353 
354  // set up the MuonResidualsFitters (which also collect residuals for fitting)
355  m_me11map.clear();
356  m_fitters.clear();
357  m_indexes.clear();
358  m_fitterOrder.clear();
359 
360  for (const auto& ali : m_alignables) {
361  bool made_fitter = false;
362 
363  // fitters for DT
364  if (ali->alignableObjectId() == align::AlignableDTChamber) {
365  DTChamberId id(ali->geomDetId().rawId());
366 
367  if (id.station() == 4) {
368  m_fitters[ali] = new MuonResidualsTwoBin(
369  m_twoBin,
370  new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment),
371  new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment));
372  made_fitter = true;
373  } else {
374  m_fitters[ali] = new MuonResidualsTwoBin(
375  m_twoBin,
376  new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment),
377  new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment));
378  made_fitter = true;
379  }
380  }
381 
382  // fitters for CSC
383  else if (ali->alignableObjectId() == align::AlignableCSCChamber) {
384  auto thisali = ali;
385  CSCDetId id(ali->geomDetId().rawId());
386 
387  // take care of ME1/1a
388  if (m_combineME11 && id.station() == 1 && id.ring() == 4) {
389  CSCDetId pairid(id.endcap(), 1, 1, id.chamber());
390 
391  for (const auto& ali2 : m_alignables) {
392  if (ali2->alignableObjectId() == align::AlignableCSCChamber && ali2->geomDetId().rawId() == pairid.rawId()) {
393  thisali = ali2;
394  break;
395  }
396  }
397  m_me11map[ali] = thisali; // points from each ME1/4 chamber to the corresponding ME1/1 chamber
398  }
399 
400  if (thisali == ali) // don't make fitters for ME1/4; they get taken care of in ME1/1
401  {
402  m_fitters[ali] = new MuonResidualsTwoBin(
403  m_twoBin,
404  new MuonResiduals6DOFrphiFitter(
405  residualsModel, m_minAlignmentHits, useResiduals, cscGeometry, m_weightAlignment),
406  new MuonResiduals6DOFrphiFitter(
407  residualsModel, m_minAlignmentHits, useResiduals, cscGeometry, m_weightAlignment));
408  made_fitter = true;
409  }
410  }
411 
412  else {
413  throw cms::Exception("MuonAlignmentFromReference")
414  << "only DTChambers and CSCChambers can be aligned with this module" << std::endl;
415  }
416 
417  if (made_fitter) {
418  m_fitters[ali]->setStrategy(m_strategy);
419 
420  int index = ali->geomDetId().rawId();
421  m_indexes.push_back(index);
422  m_fitterOrder[index] = m_fitters[ali];
423  }
424  } // end loop over chambers chosen for alignment
425 
426  // cannonical order of fitters in the file
427  std::sort(m_indexes.begin(), m_indexes.end());
428 
429  // de-weight all chambers but the reference
430  const auto& all_DT_chambers = alignableMuon->DTChambers();
431  const auto& all_CSC_chambers = alignableMuon->CSCChambers();
432  align::Alignables reference;
433  if (!m_reference.empty())
434  parseReference(reference, all_DT_chambers, all_CSC_chambers);
435 
436  alignmentParameterStore->setAlignmentPositionError(all_DT_chambers, 100000000., 0.);
437  alignmentParameterStore->setAlignmentPositionError(all_CSC_chambers, 100000000., 0.);
438  alignmentParameterStore->setAlignmentPositionError(reference, 0., 0.);
439 }
440 
441 void MuonAlignmentFromReference::run(const edm::EventSetup& iSetup, const EventInfo& eventInfo) {
442  if (m_debug)
443  std::cout << "****** EVENT START *******" << std::endl;
444  m_counter_events++;
445 
446  const GlobalTrackingGeometry* globalGeometry = &iSetup.getData(m_globTackingToken);
447  const MagneticField* magneticField = &iSetup.getData(m_MagFieldToken);
448  const Propagator* prop = &iSetup.getData(m_propToken);
449  const DetIdAssociator* muonDetIdAssociator = &iSetup.getData(m_DetIdToken);
450  auto builder = iSetup.getHandle(m_builderToken);
451 
452  if (m_muonCollectionTag.label().empty()) // use trajectories
453  {
454  if (m_debug)
455  std::cout << "JUST BEFORE LOOP OVER trajTrackPairs" << std::endl;
456  // const ConstTrajTrackPairCollection &trajtracks = eventInfo.trajTrackPairs_; // trajTrackPairs_ now private
457  const ConstTrajTrackPairCollection& trajtracks = eventInfo.trajTrackPairs();
458 
459  for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
460  ++trajtrack) {
461  m_counter_tracks++;
462 
463  const Trajectory* traj = (*trajtrack).first;
464  const reco::Track* track = (*trajtrack).second;
465 
466  if (m_minTrackPt < track->pt() && track->pt() < m_maxTrackPt && m_minTrackP < track->p() &&
467  track->p() < m_maxTrackP) {
468  m_counter_trackmomentum++;
469 
470  if (fabs(track->dxy(eventInfo.beamSpot().position())) < m_maxDxy) {
471  m_counter_trackdxy++;
472  if (m_debug)
473  std::cout << "JUST BEFORE muonResidualsFromTrack" << std::endl;
474  MuonResidualsFromTrack muonResidualsFromTrack(builder,
475  magneticField,
476  globalGeometry,
477  muonDetIdAssociator,
478  prop,
479  traj,
480  track,
481  m_alignableNavigator,
482  1000.);
483  if (m_debug)
484  std::cout << "JUST AFTER muonResidualsFromTrack" << std::endl;
485 
486  if (m_debug)
487  std::cout << "JUST BEFORE PROCESS" << std::endl;
488  processMuonResidualsFromTrack(muonResidualsFromTrack);
489  if (m_debug)
490  std::cout << "JUST AFTER PROCESS" << std::endl;
491  }
492  } // end if track p is within range
493  } // end if track pT is within range
494  if (m_debug)
495  std::cout << "JUST AFTER LOOP OVER trajTrackPairs" << std::endl;
496 
497  } else // use muons
498  {
499  /*
500  for (reco::MuonCollection::const_iterator muon = eventInfo.muonCollection_->begin(); muon != eventInfo.muonCollection_->end(); ++muon)
501  {
502  if ( !(muon->isTrackerMuon() && muon->innerTrack().isNonnull() ) ) continue;
503 
504  m_counter_tracks++;
505 
506  if (m_minTrackPt < muon->pt() && muon->pt() < m_maxTrackPt && m_minTrackP < muon->p() && muon->p() < m_maxTrackP)
507  {
508  m_counter_trackmomentum++;
509 
510  if (fabs(muon->innerTrack()->dxy(eventInfo.beamSpot_.position())) < m_maxDxy)
511  {
512  m_counter_trackdxy++;
513 
514  //std::cout<<" *** will make MuonResidualsFromTrack ***"<<std::endl;
515  MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), m_alignableNavigator, 100.);
516  //std::cout<<" *** have made MuonResidualsFromTrack ***"<<std::endl;
517 
518  //std::cout<<" trk eta="<<muon->eta()<<" ndof="<<muon->innerTrack()->ndof()<<" nchi2="<<muon->innerTrack()->normalizedChi2()
519  // <<" muresnchi2="<<muonResidualsFromTrack.normalizedChi2()<<" muresnhits="<<muonResidualsFromTrack.trackerNumHits()<<std::endl;
520 
521  processMuonResidualsFromTrack(muonResidualsFromTrack);
522  } // end if track p is within range
523  } // end if track pT is within range
524  } // end loop over tracks
525  */
526  }
527 }
528 
529 void MuonAlignmentFromReference::processMuonResidualsFromTrack(MuonResidualsFromTrack& mrft) {
530  // std::cout << "minTrackerHits: " << mrft.trackerNumHits() << std::endl;
531  if (mrft.trackerNumHits() >= m_minTrackerHits) {
532  m_counter_trackerhits++;
533  // std::cout << "mrft.normalizedChi2(): " << mrft.normalizedChi2() << std::endl;
534 
535  if (mrft.normalizedChi2() < m_maxTrackerRedChi2) {
536  m_counter_trackerchi2++;
537  if (m_allowTIDTEC || !mrft.contains_TIDTEC()) {
538  m_counter_trackertidtec++;
539 
540  std::vector<DetId> chamberIds = mrft.chamberIds();
541 
542  if ((int)chamberIds.size() >= m_minNCrossedChambers) {
543  m_counter_minchambers++;
544 
545  char charge = (mrft.getTrack()->charge() > 0 ? 1 : -1);
546 
547  for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end();
548  ++chamberId) {
549  if (chamberId->det() != DetId::Muon)
550  continue;
551  m_counter_totchambers++;
552 
553  // DT station 1,2,3
554  if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT && DTChamberId(chamberId->rawId()).station() != 4) {
555  MuonChamberResidual* dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
556  MuonChamberResidual* dt2 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT2);
557 
558  m_counter_station123++;
559  if (dt13 != nullptr && dt2 != nullptr) {
560  m_counter_station123valid++;
561  if (dt13->numHits() >= m_minDT13Hits) {
562  m_counter_station123dt13hits++;
563  if (dt2->numHits() >= m_minDT2Hits) {
564  m_counter_station123dt2hits++;
565  std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter =
566  m_fitters.find(dt13->chamberAlignable());
567  if (fitter != m_fitters.end()) {
568  m_counter_station123aligning++;
569  if (fabs(dt2->resslope()) < m_maxResSlopeY && (dt2->chi2() / double(dt2->ndof())) < 2.0) {
570  m_counter_resslopey++;
571  double* residdata = new double[MuonResiduals6DOFFitter::kNData];
572  residdata[MuonResiduals6DOFFitter::kResidX] = dt13->residual();
573  residdata[MuonResiduals6DOFFitter::kResidY] = dt2->residual();
574  residdata[MuonResiduals6DOFFitter::kResSlopeX] = dt13->resslope();
575  residdata[MuonResiduals6DOFFitter::kResSlopeY] = dt2->resslope();
576  residdata[MuonResiduals6DOFFitter::kPositionX] = dt13->trackx();
577  residdata[MuonResiduals6DOFFitter::kPositionY] = dt13->tracky();
578  residdata[MuonResiduals6DOFFitter::kAngleX] = dt13->trackdxdz();
579  residdata[MuonResiduals6DOFFitter::kAngleY] = dt13->trackdydz();
580  residdata[MuonResiduals6DOFFitter::kRedChi2] =
581  (dt13->chi2() + dt2->chi2()) / double(dt13->ndof() + dt2->ndof());
582  residdata[MuonResiduals6DOFFitter::kPz] = mrft.getTrack()->pz();
583  residdata[MuonResiduals6DOFFitter::kPt] = mrft.getTrack()->pt();
584  residdata[MuonResiduals6DOFFitter::kCharge] = mrft.getTrack()->charge();
585  residdata[MuonResiduals6DOFFitter::kStation] = DTChamberId(chamberId->rawId()).station();
586  residdata[MuonResiduals6DOFFitter::kWheel] = DTChamberId(chamberId->rawId()).wheel();
587  residdata[MuonResiduals6DOFFitter::kSector] = DTChamberId(chamberId->rawId()).sector();
588  residdata[MuonResiduals6DOFFitter::kChambW] = dt13->ChambW();
589  residdata[MuonResiduals6DOFFitter::kChambl] = dt13->Chambl();
590 
591  if (m_debug) {
592  std::cout << "processMuonResidualsFromTrack 6DOF dt13->residual() " << dt13->residual()
593  << std::endl;
594  std::cout << " dt2->residual() " << dt2->residual()
595  << std::endl;
596  std::cout << " dt13->resslope() " << dt13->resslope()
597  << std::endl;
598  std::cout << " dt2->resslope() " << dt2->resslope()
599  << std::endl;
600  std::cout << " dt13->trackx() " << dt13->trackx()
601  << std::endl;
602  std::cout << " dt13->tracky() " << dt13->tracky()
603  << std::endl;
604  std::cout << " dt13->trackdxdz() " << dt13->trackdxdz()
605  << std::endl;
606  std::cout << " dt13->trackdydz() " << dt13->trackdydz()
607  << std::endl;
608  }
609 
610  fitter->second->fill(charge, residdata);
611  // the MuonResidualsFitter will delete the array when it is destroyed
612  }
613  }
614  }
615  }
616  }
617  }
618 
619  // DT 4th station
620  else if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT &&
621  DTChamberId(chamberId->rawId()).station() == 4) {
622  MuonChamberResidual* dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
623 
624  m_counter_station4++;
625  if (dt13 != nullptr) {
626  m_counter_station4valid++;
627  if (dt13->numHits() >= m_minDT13Hits) {
628  m_counter_station4hits++;
629 
630  std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter =
631  m_fitters.find(dt13->chamberAlignable());
632  if (fitter != m_fitters.end()) {
633  m_counter_station4aligning++;
634 
635  double* residdata = new double[MuonResiduals5DOFFitter::kNData];
636  residdata[MuonResiduals5DOFFitter::kResid] = dt13->residual();
637  residdata[MuonResiduals5DOFFitter::kResSlope] = dt13->resslope();
638  residdata[MuonResiduals5DOFFitter::kPositionX] = dt13->trackx();
639  residdata[MuonResiduals5DOFFitter::kPositionY] = dt13->tracky();
640  residdata[MuonResiduals5DOFFitter::kAngleX] = dt13->trackdxdz();
641  residdata[MuonResiduals5DOFFitter::kAngleY] = dt13->trackdydz();
642  residdata[MuonResiduals5DOFFitter::kRedChi2] = dt13->chi2() / double(dt13->ndof());
643  residdata[MuonResiduals5DOFFitter::kPz] = mrft.getTrack()->pz();
644  residdata[MuonResiduals5DOFFitter::kPt] = mrft.getTrack()->pt();
645  residdata[MuonResiduals5DOFFitter::kCharge] = mrft.getTrack()->charge();
646  residdata[MuonResiduals5DOFFitter::kStation] = DTChamberId(chamberId->rawId()).station();
647  residdata[MuonResiduals5DOFFitter::kWheel] = DTChamberId(chamberId->rawId()).wheel();
648  residdata[MuonResiduals5DOFFitter::kSector] = DTChamberId(chamberId->rawId()).sector();
649  residdata[MuonResiduals5DOFFitter::kChambW] = dt13->ChambW();
650  residdata[MuonResiduals5DOFFitter::kChambl] = dt13->Chambl();
651 
652  if (m_debug) {
653  std::cout << "processMuonResidualsFromTrack 5DOF dt13->residual() " << dt13->residual()
654  << std::endl;
655  std::cout << " dt13->resslope() " << dt13->resslope()
656  << std::endl;
657  std::cout << " dt13->trackx() " << dt13->trackx()
658  << std::endl;
659  std::cout << " dt13->tracky() " << dt13->tracky()
660  << std::endl;
661  std::cout << " dt13->trackdxdz() " << dt13->trackdxdz()
662  << std::endl;
663  std::cout << " dt13->trackdydz() " << dt13->trackdydz()
664  << std::endl;
665  }
666 
667  fitter->second->fill(charge, residdata);
668  // the MuonResidualsFitter will delete the array when it is destroyed
669  }
670  }
671  }
672  } // end DT 4th station
673 
674  // CSC
675  else if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC) {
676  MuonChamberResidual* csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
677  m_counter_csc++;
678  if (csc != nullptr) {
679  m_counter_cscvalid++;
680  if (csc->numHits() >= m_minCSCHits) {
681  m_counter_cschits++;
682  Alignable* ali = csc->chamberAlignable();
683 
684  CSCDetId id(ali->geomDetId().rawId());
685  if (m_combineME11 && id.station() == 1 && id.ring() == 4)
686  ali = m_me11map[ali];
687 
688  std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(ali);
689  if (fitter != m_fitters.end()) {
690  m_counter_cscaligning++;
691  double* residdata = new double[MuonResiduals6DOFrphiFitter::kNData];
692  residdata[MuonResiduals6DOFrphiFitter::kResid] = csc->residual();
693  residdata[MuonResiduals6DOFrphiFitter::kResSlope] = csc->resslope();
694  residdata[MuonResiduals6DOFrphiFitter::kPositionX] = csc->trackx();
695  residdata[MuonResiduals6DOFrphiFitter::kPositionY] = csc->tracky();
696  residdata[MuonResiduals6DOFrphiFitter::kAngleX] = csc->trackdxdz();
697  residdata[MuonResiduals6DOFrphiFitter::kAngleY] = csc->trackdydz();
698  residdata[MuonResiduals6DOFrphiFitter::kRedChi2] = csc->chi2() / double(csc->ndof());
699  residdata[MuonResiduals6DOFrphiFitter::kPz] = mrft.getTrack()->pz();
700  residdata[MuonResiduals6DOFrphiFitter::kPt] = mrft.getTrack()->pt();
701  residdata[MuonResiduals6DOFrphiFitter::kCharge] = mrft.getTrack()->charge();
702 
703  if (m_debug) {
704  std::cout << "processMuonResidualsFromTrack 6DOFrphi csc->residual() " << csc->residual()
705  << std::endl;
706  std::cout << " csc->resslope() " << csc->resslope()
707  << std::endl;
708  std::cout << " csc->trackx() " << csc->trackx()
709  << std::endl;
710  std::cout << " csc->tracky() " << csc->tracky()
711  << std::endl;
712  std::cout << " csc->trackdxdz() " << csc->trackdxdz()
713  << std::endl;
714  std::cout << " csc->trackdydz() " << csc->trackdydz()
715  << std::endl;
716  }
717 
718  fitter->second->fill(charge, residdata);
719  // the MuonResidualsFitter will delete the array when it is destroyed
720  }
721  }
722  }
723  } // end CSC
724 
725  else if (m_doDT && m_doCSC)
726  assert(false);
727 
728  } // end loop over chamberIds
729  } // # crossed muon chambers ok
730  } // endcap tracker ok
731  } // chi2 ok
732  } // trackerNumHits ok
733 }
734 
735 void MuonAlignmentFromReference::terminate(const edm::EventSetup& iSetup) {
736  bool m_debug = false;
737 
738  // one-time print-out
739  std::cout << "Counters:" << std::endl
740  << "COUNT{ events: " << m_counter_events << " }" << std::endl
741  << "COUNT{ tracks: " << m_counter_tracks << " }" << std::endl
742  << "COUNT{ trackppt: " << m_counter_trackmomentum << " }" << std::endl
743  << "COUNT{ trackdxy: " << m_counter_trackdxy << " }" << std::endl
744  << "COUNT{ trackerhits: " << m_counter_trackerhits << " }" << std::endl
745  << "COUNT{ trackerchi2: " << m_counter_trackerchi2 << " }" << std::endl
746  << "COUNT{ trackertidtec: " << m_counter_trackertidtec << " }" << std::endl
747  << "COUNT{ minnchambers: " << m_counter_minchambers << " }" << std::endl
748  << "COUNT{ totchambers: " << m_counter_totchambers << " }" << std::endl
749  << "COUNT{ station123: " << m_counter_station123 << " }" << std::endl
750  << "COUNT{ station123valid: " << m_counter_station123valid << " }" << std::endl
751  << "COUNT{ station123dt13hits: " << m_counter_station123dt13hits << " }" << std::endl
752  << "COUNT{ station123dt2hits: " << m_counter_station123dt2hits << " }" << std::endl
753  << "COUNT{ station123aligning: " << m_counter_station123aligning << " }" << std::endl
754  << "COUNT{ resslopey: " << m_counter_resslopey << " }" << std::endl
755  << "COUNT{ station4: " << m_counter_station4 << " }" << std::endl
756  << "COUNT{ station4valid: " << m_counter_station4valid << " }" << std::endl
757  << "COUNT{ station4hits: " << m_counter_station4hits << " }" << std::endl
758  << "COUNT{ station4aligning: " << m_counter_station4aligning << " }" << std::endl
759  << "COUNT{ csc: " << m_counter_csc << " }" << std::endl
760  << "COUNT{ cscvalid: " << m_counter_cscvalid << " }" << std::endl
761  << "COUNT{ cschits: " << m_counter_cschits << " }" << std::endl
762  << "COUNT{ cscaligning: " << m_counter_cscaligning << " }" << std::endl
763  << "That's all!" << std::endl;
764 
765  TStopwatch stop_watch;
766 
767  // collect temporary files
768  if (!m_readTemporaryFiles.empty()) {
769  stop_watch.Start();
770  readTmpFiles();
771  if (m_debug)
772  std::cout << "readTmpFiles took " << stop_watch.CpuTime() << " sec" << std::endl;
773  stop_watch.Stop();
774  }
775 
776  // select residuals peaks and discard tails if peakNSigma>0 (only while doing alignment)
777  if (m_peakNSigma > 0. && m_doAlignment) {
778  stop_watch.Start();
779  selectResidualsPeaks();
780  if (m_debug)
781  std::cout << "selectResidualsPeaks took " << stop_watch.CpuTime() << " sec" << std::endl;
782  stop_watch.Stop();
783  }
784 
785  if (m_BFieldCorrection > 0 && m_doAlignment) {
786  stop_watch.Start();
787  correctBField();
788  if (m_debug)
789  std::cout << "correctBField took " << stop_watch.CpuTime() << " sec" << std::endl;
790  stop_watch.Stop();
791  }
792 
793  if (m_doAlignment && !m_doCSC) // for now apply fiducial cuts to DT only
794  {
795  stop_watch.Start();
796  fiducialCuts();
797  if (m_debug)
798  std::cout << "fiducialCuts took " << stop_watch.CpuTime() << " sec" << std::endl;
799  stop_watch.Stop();
800  }
801 
802  // optionally, create an nutuple for easy debugging
803  if (m_createNtuple) {
804  stop_watch.Start();
805  fillNtuple();
806  if (m_debug)
807  std::cout << "fillNtuple took " << stop_watch.CpuTime() << " sec" << std::endl;
808  stop_watch.Stop();
809  }
810 
811  if (m_doAlignment) {
812  stop_watch.Start();
813  eraseNotSelectedResiduals();
814  if (m_debug)
815  std::cout << "eraseNotSelectedResiduals took " << stop_watch.CpuTime() << " sec" << std::endl;
816  stop_watch.Stop();
817  }
818 
819  // fit and align (time-consuming, so the user can turn it off if in a residuals-gathering job)
820  if (m_doAlignment) {
821  stop_watch.Start();
822  fitAndAlign();
823  if (m_debug)
824  std::cout << "fitAndAlign took " << stop_watch.CpuTime() << " sec" << std::endl;
825  stop_watch.Stop();
826  }
827 
828  // write out the pseudontuples for a later job to collect
829  if (m_writeTemporaryFile != std::string(""))
830  writeTmpFiles();
831  if (m_debug)
832  std::cout << "end: MuonAlignmentFromReference::terminate()" << std::endl;
833 }
834 
835 void MuonAlignmentFromReference::fitAndAlign() {
836  bool m_debug = false;
837 
838  edm::Service<TFileService> tfileService;
839  TFileDirectory rootDirectory(tfileService->mkdir("MuonAlignmentFromReference"));
840 
841  std::ofstream report;
842  bool writeReport = (m_reportFileName != std::string(""));
843  if (writeReport) {
844  report.open(m_reportFileName.c_str());
845  report << "nan = None; NAN = None" << std::endl;
846  report << "nan = 0" << std::endl;
847  report << "reports = []" << std::endl;
848  report << "class ValErr:" << std::endl
849  << " def __init__(self, value, error, antisym):" << std::endl
850  << " self.value, self.error, self.antisym = value, error, antisym" << std::endl
851  << "" << std::endl
852  << " def __repr__(self):" << std::endl
853  << " if self.antisym == 0.:" << std::endl
854  << " return \"%g +- %g\" % (self.value, self.error)" << std::endl
855  << " else:" << std::endl
856  << " return \"%g +- %g ~ %g\" % (self.value, self.error, self.antisym)" << std::endl
857  << "" << std::endl
858  << "class Report:" << std::endl
859  << " def __init__(self, chamberId, postal_address, name):" << std::endl
860  << " self.chamberId, self.postal_address, self.name = chamberId, postal_address, name" << std::endl
861  << " self.status = \"NOFIT\"" << std::endl
862  << " self.fittype = None" << std::endl
863  << "" << std::endl
864  << " def add_parameters(self, deltax, deltay, deltaz, deltaphix, deltaphiy, deltaphiz, loglikelihood, "
865  "numsegments, sumofweights, redchi2):"
866  << std::endl
867  << " self.status = \"PASS\"" << std::endl
868  << " self.deltax, self.deltay, self.deltaz, self.deltaphix, self.deltaphiy, self.deltaphiz = deltax, "
869  "deltay, deltaz, deltaphix, deltaphiy, deltaphiz"
870  << std::endl
871  << " self.loglikelihood, self.numsegments, self.sumofweights, self.redchi2 = loglikelihood, "
872  "numsegments, sumofweights, redchi2"
873  << std::endl
874  << "" << std::endl
875  << " def add_stats(self, median_x, median_y, median_dxdz, median_dydz, mean30_x, mean30_y, mean20_dxdz, "
876  "mean50_dydz, mean15_x, mean15_y, mean10_dxdz, mean25_dydz, wmean30_x, wmean30_y, wmean20_dxdz, "
877  "wmean50_dydz, wmean15_x, wmean15_y, wmean10_dxdz, wmean25_dydz, stdev30_x, stdev30_y, stdev20_dxdz, "
878  "stdev50_dydz, stdev15_x, stdev15_y, stdev10_dxdz, stdev25_dydz):"
879  << std::endl
880  << " self.median_x, self.median_y, self.median_dxdz, self.median_dydz, self.mean30_x, self.mean30_y, "
881  "self.mean20_dxdz, self.mean50_dydz, self.mean15_x, self.mean15_y, self.mean10_dxdz, self.mean25_dydz, "
882  "self.wmean30_x, self.wmean30_y, self.wmean20_dxdz, self.wmean50_dydz, self.wmean15_x, self.wmean15_y, "
883  "self.wmean10_dxdz, self.wmean25_dydz, self.stdev30_x, self.stdev30_y, self.stdev20_dxdz, "
884  "self.stdev50_dydz, self.stdev15_x, self.stdev15_y, self.stdev10_dxdz, self.stdev25_dydz = median_x, "
885  "median_y, median_dxdz, median_dydz, mean30_x, mean30_y, mean20_dxdz, mean50_dydz, mean15_x, mean15_y, "
886  "mean10_dxdz, mean25_dydz, wmean30_x, wmean30_y, wmean20_dxdz, wmean50_dydz, wmean15_x, wmean15_y, "
887  "wmean10_dxdz, wmean25_dydz, stdev30_x, stdev30_y, stdev20_dxdz, stdev50_dydz, stdev15_x, stdev15_y, "
888  "stdev10_dxdz, stdev25_dydz"
889  << std::endl
890  << "" << std::endl
891  << " def __repr__(self):" << std::endl
892  << " return \"<Report %s %s %s>\" % (self.postal_address[0], \" \".join(map(str, "
893  "self.postal_address[1:])), self.status)"
894  << std::endl
895  << std::endl;
896  }
897 
898  if (m_debug)
899  std::cout << "***** just after report.open" << std::endl;
900 
901  for (const auto& ali : m_alignables) {
902  if (m_debug)
903  std::cout << "***** Start loop over alignables" << std::endl;
904 
905  std::vector<bool> selector = ali->alignmentParameters()->selector();
906  bool align_x = selector[0];
907  bool align_y = selector[1];
908  bool align_z = selector[2];
909  bool align_phix = selector[3];
910  bool align_phiy = selector[4];
911  bool align_phiz = selector[5];
912  int numParams = ((align_x ? 1 : 0) + (align_y ? 1 : 0) + (align_z ? 1 : 0) + (align_phix ? 1 : 0) +
913  (align_phiy ? 1 : 0) + (align_phiz ? 1 : 0));
914 
915  // map from 0-5 to the index of params, above
916  std::vector<int> paramIndex;
917  int paramIndex_counter = -1;
918  if (align_x)
919  paramIndex_counter++;
920  paramIndex.push_back(paramIndex_counter);
921  if (align_y)
922  paramIndex_counter++;
923  paramIndex.push_back(paramIndex_counter);
924  if (align_z)
925  paramIndex_counter++;
926  paramIndex.push_back(paramIndex_counter);
927  if (align_phix)
928  paramIndex_counter++;
929  paramIndex.push_back(paramIndex_counter);
930  if (align_phiy)
931  paramIndex_counter++;
932  paramIndex.push_back(paramIndex_counter);
933  if (align_phiz)
934  paramIndex_counter++;
935  paramIndex.push_back(paramIndex_counter);
936 
937  DetId id = ali->geomDetId();
938 
939  auto thisali = ali;
940  if (m_combineME11 && id.subdetId() == MuonSubdetId::CSC) {
941  CSCDetId cscid(id.rawId());
942  if (cscid.station() == 1 && cscid.ring() == 4)
943  thisali = m_me11map[ali];
944  }
945 
946  if (m_debug)
947  std::cout << "***** loop over alignables 1" << std::endl;
948 
949  char cname[40];
950  char wheel_label[][2] = {"A", "B", "C", "D", "E"};
951 
952  if (id.subdetId() == MuonSubdetId::DT) {
953  DTChamberId chamberId(id.rawId());
954 
955  //if ( ! ( (chamberId.station()==1&&chamberId.wheel()==0) || (chamberId.station()==4&&chamberId.wheel()==2) ) ) continue;
956 
957  sprintf(cname, "MBwh%sst%dsec%02d", wheel_label[chamberId.wheel() + 2], chamberId.station(), chamberId.sector());
958  if (writeReport) {
959  report << "reports.append(Report(" << id.rawId() << ", (\"DT\", " << chamberId.wheel() << ", "
960  << chamberId.station() << ", " << chamberId.sector() << "), \"" << cname << "\"))" << std::endl;
961  }
962  } else if (id.subdetId() == MuonSubdetId::CSC) {
963  CSCDetId chamberId(id.rawId());
964  sprintf(cname,
965  "ME%s%d%d_%02d",
966  (chamberId.endcap() == 1 ? "p" : "m"),
967  chamberId.station(),
968  chamberId.ring(),
969  chamberId.chamber());
970 
971  //if ( chamberId.chamber()>6 || chamberId.endcap()==2 || ! ( (chamberId.station()==2&&chamberId.ring()==1) || (chamberId.station()==3&&chamberId.ring()==2) ) ) continue;
972 
973  if (writeReport) {
974  report << "reports.append(Report(" << id.rawId() << ", (\"CSC\", " << chamberId.endcap() << ", "
975  << chamberId.station() << ", " << chamberId.ring() << ", " << chamberId.chamber() << "), \"" << cname
976  << "\"))" << std::endl;
977  }
978  }
979 
980  if (m_debug)
981  std::cout << "***** loop over alignables 2" << std::endl;
982 
983  //if(! ( strcmp(cname,"MBwhCst3sec12")==0 || strcmp(cname,"MBwhCst3sec06")==0)) continue;
984 
985  std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(thisali);
986 
987  if (m_debug)
988  std::cout << "***** loop over alignables 3" << std::endl;
989 
990  if (fitter != m_fitters.end()) {
991  //if (fitter->second->type() != MuonResidualsFitter::k6DOFrphi) continue;
992 
993  TStopwatch stop_watch;
994  stop_watch.Start();
995 
996  // MINUIT is verbose in std::cout anyway
997  if (m_debug)
998  std::cout << "============================================================================================="
999  << std::endl;
1000  if (m_debug)
1001  std::cout << "Fitting " << cname << std::endl;
1002 
1003  if (writeReport) {
1004  report << "reports[-1].posNum = " << fitter->second->numResidualsPos() << std::endl;
1005  report << "reports[-1].negNum = " << fitter->second->numResidualsNeg() << std::endl;
1006  }
1007 
1008  if (fitter->second->type() == MuonResidualsFitter::k5DOF) {
1009  if (!align_x)
1010  fitter->second->fix(MuonResiduals5DOFFitter::kAlignX);
1011  if (!align_z)
1012  fitter->second->fix(MuonResiduals5DOFFitter::kAlignZ);
1013  if (!align_phix)
1014  fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiX);
1015  if (!align_phiy)
1016  fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiY);
1017  if (!align_phiz)
1018  fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiZ);
1019  } else if (fitter->second->type() == MuonResidualsFitter::k6DOF) {
1020  if (!align_x)
1021  fitter->second->fix(MuonResiduals6DOFFitter::kAlignX);
1022  if (!align_y)
1023  fitter->second->fix(MuonResiduals6DOFFitter::kAlignY);
1024  if (!align_z)
1025  fitter->second->fix(MuonResiduals6DOFFitter::kAlignZ);
1026  if (!align_phix)
1027  fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiX);
1028  if (!align_phiy)
1029  fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiY);
1030  if (!align_phiz)
1031  fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiZ);
1032  } else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi) {
1033  if (!align_x)
1034  fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignX);
1035  if (!align_y)
1036  fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignY);
1037  if (!align_z)
1038  fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignZ);
1039  if (!align_phix)
1040  fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1041  if (!align_phiy)
1042  fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1043  if (!align_phiz)
1044  fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1045  } else
1046  assert(false);
1047 
1048  if (m_debug)
1049  std::cout << "***** loop over alignables 4" << std::endl;
1050 
1051  AlgebraicVector params(numParams);
1052  AlgebraicSymMatrix cov(numParams);
1053 
1054  if (fitter->second->numsegments() >= m_minAlignmentHits) {
1055  if (m_debug)
1056  std::cout << "***** loop over alignables 5" << std::endl;
1057 
1058  bool successful_fit = fitter->second->fit(thisali);
1059 
1060  if (m_debug)
1061  std::cout << "***** loop over alignables 6 " << fitter->second->type() << std::endl;
1062 
1063  double loglikelihood = fitter->second->loglikelihood();
1064  double numsegments = fitter->second->numsegments();
1065  double sumofweights = fitter->second->sumofweights();
1066  double redchi2 = fitter->second->plot(cname, &rootDirectory, thisali);
1067 
1068  if (fitter->second->type() == MuonResidualsFitter::k5DOF) {
1069  if (m_debug)
1070  std::cout << "***** loop over alignables k5DOF" << std::endl;
1071 
1072  double deltax_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignX);
1073  double deltax_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignX);
1074  double deltax_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignX);
1075 
1076  double deltaz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignZ);
1077  double deltaz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignZ);
1078  double deltaz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignZ);
1079 
1080  double deltaphix_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiX);
1081  double deltaphix_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiX);
1082  double deltaphix_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiX);
1083 
1084  double deltaphiy_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiY);
1085  double deltaphiy_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiY);
1086  double deltaphiy_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiY);
1087 
1088  double deltaphiz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiZ);
1089  double deltaphiz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiZ);
1090  double deltaphiz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiZ);
1091 
1092  double sigmaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidSigma);
1093  double sigmaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidSigma);
1094  double sigmaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidSigma);
1095 
1096  double sigmaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeSigma);
1097  double sigmaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeSigma);
1098  double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeSigma);
1099 
1100  double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error,
1101  gammaresslope_antisym;
1102  gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error =
1103  gammaresslope_antisym = 0.;
1104 
1105  if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1106  fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1107  fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1108  gammaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidGamma);
1109  gammaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidGamma);
1110  gammaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidGamma);
1111 
1112  gammaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeGamma);
1113  gammaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeGamma);
1114  gammaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeGamma);
1115  }
1116 
1117  if (writeReport) {
1118  report << "reports[-1].fittype = \"5DOF\"" << std::endl;
1119  report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", "
1120  << deltax_antisym << "), \\" << std::endl
1121  << " None, \\" << std::endl
1122  << " ValErr(" << deltaz_value << ", " << deltaz_error << ", "
1123  << deltaz_antisym << "), \\" << std::endl
1124  << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", "
1125  << deltaphix_antisym << "), \\" << std::endl
1126  << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", "
1127  << deltaphiy_antisym << "), \\" << std::endl
1128  << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", "
1129  << deltaphiz_antisym << "), \\" << std::endl
1130  << " " << loglikelihood << ", " << numsegments << ", " << sumofweights
1131  << ", " << redchi2 << ")" << std::endl;
1132  report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", "
1133  << sigmaresid_antisym << ")" << std::endl;
1134  report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error
1135  << ", " << sigmaresslope_antisym << ")" << std::endl;
1136  if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1137  fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1138  fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1139  report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", "
1140  << gammaresid_antisym << ")" << std::endl;
1141  report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error
1142  << ", " << gammaresslope_antisym << ")" << std::endl;
1143  }
1144 
1145  report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals5DOFFitter::kResid) << ", "
1146  << "None, " << fitter->second->median(MuonResiduals5DOFFitter::kResSlope) << ", "
1147  << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 30.) << ", "
1148  << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 20.) << ", "
1149  << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 15.) << ", "
1150  << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 10.) << ", "
1151  << "None, "
1152  << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 30.)
1153  << ", "
1154  << "None, "
1155  << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 20.)
1156  << ", "
1157  << "None, "
1158  << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 15.)
1159  << ", "
1160  << "None, "
1161  << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 10.)
1162  << ", "
1163  << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 30.) << ", "
1164  << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 20.) << ", "
1165  << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 15.) << ", "
1166  << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 10.) << ", "
1167  << "None)" << std::endl;
1168 
1169  std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
1170  namesimple_x << cname << "_simple_x";
1171  namesimple_dxdz << cname << "_simple_dxdz";
1172  nameweighted_x << cname << "_weighted_x";
1173  nameweighted_dxdz << cname << "_weighted_dxdz";
1174 
1175  fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals5DOFFitter::kResid, 10.);
1176  fitter->second->plotsimple(
1177  namesimple_dxdz.str(), &rootDirectory, MuonResiduals5DOFFitter::kResSlope, 1000.);
1178 
1179  fitter->second->plotweighted(nameweighted_x.str(),
1180  &rootDirectory,
1181  MuonResiduals5DOFFitter::kResid,
1182  MuonResiduals5DOFFitter::kRedChi2,
1183  10.);
1184  fitter->second->plotweighted(nameweighted_dxdz.str(),
1185  &rootDirectory,
1186  MuonResiduals5DOFFitter::kResSlope,
1187  MuonResiduals5DOFFitter::kRedChi2,
1188  1000.);
1189  }
1190 
1191  if (successful_fit) {
1192  if (align_x)
1193  params[paramIndex[0]] = deltax_value;
1194  if (align_z)
1195  params[paramIndex[2]] = deltaz_value;
1196  if (align_phix)
1197  params[paramIndex[3]] = deltaphix_value;
1198  if (align_phiy)
1199  params[paramIndex[4]] = deltaphiy_value;
1200  if (align_phiz)
1201  params[paramIndex[5]] = deltaphiz_value;
1202  }
1203  } // end if 5DOF
1204 
1205  else if (fitter->second->type() == MuonResidualsFitter::k6DOF) {
1206  if (m_debug)
1207  std::cout << "***** loop over alignables k6DOF" << std::endl;
1208 
1209  double deltax_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignX);
1210  double deltax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignX);
1211  double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignX);
1212 
1213  double deltay_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignY);
1214  double deltay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignY);
1215  double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignY);
1216 
1217  double deltaz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignZ);
1218  double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignZ);
1219  double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignZ);
1220 
1221  double deltaphix_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiX);
1222  double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiX);
1223  double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiX);
1224 
1225  double deltaphiy_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiY);
1226  double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiY);
1227  double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiY);
1228 
1229  double deltaphiz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiZ);
1230  double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiZ);
1231  double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiZ);
1232 
1233  double sigmax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXSigma);
1234  double sigmax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXSigma);
1235  double sigmax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXSigma);
1236 
1237  double sigmay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYSigma);
1238  double sigmay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYSigma);
1239  double sigmay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYSigma);
1240 
1241  double sigmadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXSigma);
1242  double sigmadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXSigma);
1243  double sigmadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXSigma);
1244 
1245  double sigmadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYSigma);
1246  double sigmadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYSigma);
1247  double sigmadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYSigma);
1248 
1249  double gammax_value, gammax_error, gammax_antisym, gammay_value, gammay_error, gammay_antisym,
1250  gammadxdz_value, gammadxdz_error, gammadxdz_antisym, gammadydz_value, gammadydz_error, gammadydz_antisym;
1251  gammax_value = gammax_error = gammax_antisym = gammay_value = gammay_error = gammay_antisym =
1252  gammadxdz_value = gammadxdz_error = gammadxdz_antisym = gammadydz_value = gammadydz_error =
1253  gammadydz_antisym = 0.;
1254  if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1255  fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1256  fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1257  gammax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXGamma);
1258  gammax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXGamma);
1259  gammax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXGamma);
1260 
1261  gammay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYGamma);
1262  gammay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYGamma);
1263  gammay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYGamma);
1264 
1265  gammadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXGamma);
1266  gammadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXGamma);
1267  gammadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXGamma);
1268 
1269  gammadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYGamma);
1270  gammadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYGamma);
1271  gammadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYGamma);
1272  }
1273 
1274  if (writeReport) {
1275  report << "reports[-1].fittype = \"6DOF\"" << std::endl;
1276  report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", "
1277  << deltax_antisym << "), \\" << std::endl
1278  << " ValErr(" << deltay_value << ", " << deltay_error << ", "
1279  << deltay_antisym << "), \\" << std::endl
1280  << " ValErr(" << deltaz_value << ", " << deltaz_error << ", "
1281  << deltaz_antisym << "), \\" << std::endl
1282  << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", "
1283  << deltaphix_antisym << "), \\" << std::endl
1284  << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", "
1285  << deltaphiy_antisym << "), \\" << std::endl
1286  << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", "
1287  << deltaphiz_antisym << "), \\" << std::endl
1288  << " " << loglikelihood << ", " << numsegments << ", " << sumofweights
1289  << ", " << redchi2 << ")" << std::endl;
1290  report << "reports[-1].sigmax = ValErr(" << sigmax_value << ", " << sigmax_error << ", " << sigmax_antisym
1291  << ")" << std::endl;
1292  report << "reports[-1].sigmay = ValErr(" << sigmay_value << ", " << sigmay_error << ", " << sigmay_antisym
1293  << ")" << std::endl;
1294  report << "reports[-1].sigmadxdz = ValErr(" << sigmadxdz_value << ", " << sigmadxdz_error << ", "
1295  << sigmadxdz_antisym << ")" << std::endl;
1296  report << "reports[-1].sigmadydz = ValErr(" << sigmadydz_value << ", " << sigmadydz_error << ", "
1297  << sigmadydz_antisym << ")" << std::endl;
1298  if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1299  fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1300  fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1301  report << "reports[-1].gammax = ValErr(" << gammax_value << ", " << gammax_error << ", " << gammax_antisym
1302  << ")" << std::endl;
1303  report << "reports[-1].gammay = ValErr(" << gammay_value << ", " << gammay_error << ", " << gammay_antisym
1304  << ")" << std::endl;
1305  report << "reports[-1].gammadxdz = ValErr(" << gammadxdz_value << ", " << gammadxdz_error << ", "
1306  << gammadxdz_antisym << ")" << std::endl;
1307  report << "reports[-1].gammadydz = ValErr(" << gammadydz_value << ", " << gammadydz_error << ", "
1308  << gammadydz_antisym << ")" << std::endl;
1309  }
1310 
1311  report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals6DOFFitter::kResidX) << ", "
1312  << fitter->second->median(MuonResiduals6DOFFitter::kResidY) << ", "
1313  << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeX) << ", "
1314  << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeY) << ", "
1315  << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
1316  << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
1317  << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
1318  << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
1319  << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
1320  << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
1321  << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
1322  << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ", "
1323  << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 30.)
1324  << ", "
1325  << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 30.)
1326  << ", "
1327  << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 20.)
1328  << ", "
1329  << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 50.)
1330  << ", "
1331  << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 15.)
1332  << ", "
1333  << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 15.)
1334  << ", "
1335  << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 10.)
1336  << ", "
1337  << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 25.)
1338  << ", " << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
1339  << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
1340  << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
1341  << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
1342  << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
1343  << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
1344  << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
1345  << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ")" << std::endl;
1346 
1347  std::stringstream namesimple_x, namesimple_y, namesimple_dxdz, namesimple_dydz, nameweighted_x,
1348  nameweighted_y, nameweighted_dxdz, nameweighted_dydz;
1349  namesimple_x << cname << "_simple_x";
1350  namesimple_y << cname << "_simple_y";
1351  namesimple_dxdz << cname << "_simple_dxdz";
1352  namesimple_dydz << cname << "_simple_dydz";
1353  nameweighted_x << cname << "_weighted_x";
1354  nameweighted_y << cname << "_weighted_y";
1355  nameweighted_dxdz << cname << "_weighted_dxdz";
1356  nameweighted_dydz << cname << "_weighted_dydz";
1357 
1358  fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidX, 10.);
1359  fitter->second->plotsimple(namesimple_y.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidY, 10.);
1360  fitter->second->plotsimple(
1361  namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeX, 1000.);
1362  fitter->second->plotsimple(
1363  namesimple_dydz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeY, 1000.);
1364 
1365  fitter->second->plotweighted(nameweighted_x.str(),
1366  &rootDirectory,
1367  MuonResiduals6DOFFitter::kResidX,
1368  MuonResiduals6DOFFitter::kRedChi2,
1369  10.);
1370  fitter->second->plotweighted(nameweighted_y.str(),
1371  &rootDirectory,
1372  MuonResiduals6DOFFitter::kResidY,
1373  MuonResiduals6DOFFitter::kRedChi2,
1374  10.);
1375  fitter->second->plotweighted(nameweighted_dxdz.str(),
1376  &rootDirectory,
1377  MuonResiduals6DOFFitter::kResSlopeX,
1378  MuonResiduals6DOFFitter::kRedChi2,
1379  1000.);
1380  fitter->second->plotweighted(nameweighted_dydz.str(),
1381  &rootDirectory,
1382  MuonResiduals6DOFFitter::kResSlopeY,
1383  MuonResiduals6DOFFitter::kRedChi2,
1384  1000.);
1385  }
1386 
1387  if (successful_fit) {
1388  if (align_x)
1389  params[paramIndex[0]] = deltax_value;
1390  if (align_y)
1391  params[paramIndex[1]] = deltay_value;
1392  if (align_z)
1393  params[paramIndex[2]] = deltaz_value;
1394  if (align_phix)
1395  params[paramIndex[3]] = deltaphix_value;
1396  if (align_phiy)
1397  params[paramIndex[4]] = deltaphiy_value;
1398  if (align_phiz)
1399  params[paramIndex[5]] = deltaphiz_value;
1400  }
1401  } // end if 6DOF
1402 
1403  else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi) {
1404  if (m_debug)
1405  std::cout << "***** loop over alignables k6DOFrphi" << std::endl;
1406 
1407  double deltax_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignX);
1408  double deltax_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignX);
1409  double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignX);
1410 
1411  double deltay_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignY);
1412  double deltay_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignY);
1413  double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignY);
1414 
1415  double deltaz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignZ);
1416  double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignZ);
1417  double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignZ);
1418 
1419  double deltaphix_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1420  double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1421  double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1422 
1423  double deltaphiy_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1424  double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1425  double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1426 
1427  double deltaphiz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1428  double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1429  double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1430 
1431  double sigmaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidSigma);
1432  double sigmaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidSigma);
1433  double sigmaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidSigma);
1434 
1435  double sigmaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
1436  double sigmaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
1437  double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
1438 
1439  double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error,
1440  gammaresslope_antisym;
1441  gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error =
1442  gammaresslope_antisym = 0.;
1443  if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1444  fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1445  fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1446  gammaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidGamma);
1447  gammaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidGamma);
1448  gammaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidGamma);
1449 
1450  gammaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
1451  gammaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
1452  gammaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
1453  }
1454 
1455  if (writeReport) {
1456  report << "reports[-1].fittype = \"6DOFrphi\"" << std::endl;
1457  report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", "
1458  << deltax_antisym << "), \\" << std::endl
1459  << " ValErr(" << deltay_value << ", " << deltay_error << ", "
1460  << deltay_antisym << "), \\" << std::endl
1461  << " ValErr(" << deltaz_value << ", " << deltaz_error << ", "
1462  << deltaz_antisym << "), \\" << std::endl
1463  << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", "
1464  << deltaphix_antisym << "), \\" << std::endl
1465  << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", "
1466  << deltaphiy_antisym << "), \\" << std::endl
1467  << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", "
1468  << deltaphiz_antisym << "), \\" << std::endl
1469  << " " << loglikelihood << ", " << numsegments << ", " << sumofweights
1470  << ", " << redchi2 << ")" << std::endl;
1471  report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", "
1472  << sigmaresid_antisym << ")" << std::endl;
1473  report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error
1474  << ", " << sigmaresslope_antisym << ")" << std::endl;
1475  if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1476  fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1477  fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1478  report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", "
1479  << gammaresid_antisym << ")" << std::endl;
1480  report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error
1481  << ", " << gammaresslope_antisym << ")" << std::endl;
1482  }
1483 
1484  report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals6DOFrphiFitter::kResid) << ", "
1485  << "None, " << fitter->second->median(MuonResiduals6DOFrphiFitter::kResSlope) << ", "
1486  << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", "
1487  << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", "
1488  << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", "
1489  << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", "
1490  << "None, "
1491  << fitter->second->wmean(
1492  MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 30.)
1493  << ", "
1494  << "None, "
1495  << fitter->second->wmean(
1496  MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 20.)
1497  << ", "
1498  << "None, "
1499  << fitter->second->wmean(
1500  MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 15.)
1501  << ", "
1502  << "None, "
1503  << fitter->second->wmean(
1504  MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 10.)
1505  << ", "
1506  << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", "
1507  << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", "
1508  << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", "
1509  << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", "
1510  << "None)" << std::endl;
1511 
1512  std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
1513  namesimple_x << cname << "_simple_x";
1514  namesimple_dxdz << cname << "_simple_dxdz";
1515  nameweighted_x << cname << "_weighted_x";
1516  nameweighted_dxdz << cname << "_weighted_dxdz";
1517 
1518  fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResid, 10.);
1519  fitter->second->plotsimple(
1520  namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResSlope, 1000.);
1521 
1522  fitter->second->plotweighted(nameweighted_x.str(),
1523  &rootDirectory,
1524  MuonResiduals6DOFrphiFitter::kResid,
1525  MuonResiduals6DOFrphiFitter::kRedChi2,
1526  10.);
1527  fitter->second->plotweighted(nameweighted_dxdz.str(),
1528  &rootDirectory,
1529  MuonResiduals6DOFrphiFitter::kResSlope,
1530  MuonResiduals6DOFrphiFitter::kRedChi2,
1531  1000.);
1532  }
1533 
1534  if (successful_fit) {
1535  if (align_x)
1536  params[paramIndex[0]] = deltax_value;
1537  if (align_y)
1538  params[paramIndex[1]] = deltay_value;
1539  if (align_z)
1540  params[paramIndex[2]] = deltaz_value;
1541  if (align_phix)
1542  params[paramIndex[3]] = deltaphix_value;
1543  if (align_phiy)
1544  params[paramIndex[4]] = deltaphiy_value;
1545  if (align_phiz)
1546  params[paramIndex[5]] = deltaphiz_value;
1547  }
1548  } // end if 6DOFrphi
1549 
1550  if (successful_fit) {
1551  align::Alignables oneortwo;
1552  oneortwo.push_back(ali);
1553  if (thisali != ali)
1554  oneortwo.push_back(thisali);
1555  m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 0., 0.);
1556  } else {
1557  if (m_debug)
1558  std::cout << "MINUIT fit failed!" << std::endl;
1559  if (writeReport) {
1560  report << "reports[-1].status = \"MINUITFAIL\"" << std::endl;
1561  }
1562 
1563  for (int i = 0; i < numParams; i++)
1564  cov[i][i] = 1000.;
1565 
1566  align::Alignables oneortwo;
1567  oneortwo.push_back(ali);
1568  if (thisali != ali)
1569  oneortwo.push_back(thisali);
1570  m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
1571  }
1572  } else { // too few hits
1573  if (m_debug)
1574  std::cout << "Too few hits!" << std::endl;
1575  if (writeReport) {
1576  report << "reports[-1].status = \"TOOFEWHITS\"" << std::endl;
1577  }
1578 
1579  for (int i = 0; i < numParams; i++)
1580  cov[i][i] = 1000.;
1581 
1582  align::Alignables oneortwo;
1583  oneortwo.push_back(ali);
1584  if (thisali != ali)
1585  oneortwo.push_back(thisali);
1586  m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
1587  }
1588 
1589  AlignmentParameters* parnew = ali->alignmentParameters()->cloneFromSelected(params, cov);
1590  ali->setAlignmentParameters(parnew);
1591  m_alignmentParameterStore->applyParameters(ali);
1592  ali->alignmentParameters()->setValid(true);
1593 
1594  if (m_debug)
1595  std::cout << cname << " fittime= " << stop_watch.CpuTime() << " sec" << std::endl;
1596  } // end we have a fitter for this alignable
1597 
1598  if (writeReport)
1599  report << std::endl;
1600 
1601  } // end loop over alignables
1602 
1603  if (writeReport)
1604  report.close();
1605 }
1606 
1607 void MuonAlignmentFromReference::readTmpFiles() {
1608  for (std::vector<std::string>::const_iterator fileName = m_readTemporaryFiles.begin();
1609  fileName != m_readTemporaryFiles.end();
1610  ++fileName) {
1611  FILE* file;
1612  int size;
1613  file = fopen((*fileName).c_str(), "r");
1614  if (file == nullptr)
1615  throw cms::Exception("MuonAlignmentFromReference")
1616  << "file \"" << *fileName << "\" can't be opened (doesn't exist?)" << std::endl;
1617 
1618  fread(&size, sizeof(int), 1, file);
1619  if (int(m_indexes.size()) != size)
1620  throw cms::Exception("MuonAlignmentFromReference")
1621  << "file \"" << *fileName << "\" has " << size << " fitters, but this job has " << m_indexes.size()
1622  << " fitters (probably corresponds to the wrong alignment job)" << std::endl;
1623 
1624  int i = 0;
1625  for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i) {
1626  MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1627  unsigned int index_toread;
1628  fread(&index_toread, sizeof(unsigned int), 1, file);
1629  if (*index != index_toread)
1630  throw cms::Exception("MuonAlignmentFromReference")
1631  << "file \"" << *fileName << "\" has index " << index_toread << " at position " << i
1632  << ", but this job is expecting " << *index << " (probably corresponds to the wrong alignment job)"
1633  << std::endl;
1634  fitter->read(file, i);
1635  }
1636 
1637  fclose(file);
1638  }
1639 }
1640 
1641 void MuonAlignmentFromReference::writeTmpFiles() {
1642  FILE* file;
1643  file = fopen(m_writeTemporaryFile.c_str(), "w");
1644  int size = m_indexes.size();
1645  fwrite(&size, sizeof(int), 1, file);
1646 
1647  int i = 0;
1648  for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i) {
1649  MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1650  unsigned int index_towrite = *index;
1651  fwrite(&index_towrite, sizeof(unsigned int), 1, file);
1652  fitter->write(file, i);
1653  }
1654 
1655  fclose(file);
1656 }
1657 
1658 void MuonAlignmentFromReference::correctBField() {
1659  bool m_debug = false;
1660 
1661  for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1662  if (m_debug)
1663  std::cout << "correcting B in " << chamberPrettyNameFromId(*index) << std::endl;
1664  MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1665  fitter->correctBField();
1666  }
1667 }
1668 
1669 void MuonAlignmentFromReference::fiducialCuts() {
1670  for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1671  if (m_debug)
1672  std::cout << "applying fiducial cuts in " << chamberPrettyNameFromId(*index) << std::endl;
1673  MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1674  fitter->fiducialCuts();
1675  }
1676 }
1677 
1678 void MuonAlignmentFromReference::eraseNotSelectedResiduals() {
1679  for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1680  if (m_debug)
1681  std::cout << "erasing in " << chamberPrettyNameFromId(*index) << std::endl;
1682  MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1683  fitter->eraseNotSelectedResiduals();
1684  }
1685 }
1686 
1687 void MuonAlignmentFromReference::selectResidualsPeaks() {
1688  // should not be called with negative peakNSigma
1689  assert(m_peakNSigma > 0.);
1690 
1691  for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1692  MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1693 
1694  int nvar = 2;
1695  int vars_index[10] = {0, 1};
1696  if (fitter->type() == MuonResidualsFitter::k5DOF) {
1697  if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 ||
1698  fitter->useRes() == MuonResidualsFitter::k1010) {
1699  nvar = 2;
1700  vars_index[0] = MuonResiduals5DOFFitter::kResid;
1701  vars_index[1] = MuonResiduals5DOFFitter::kResSlope;
1702  } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
1703  nvar = 1;
1704  vars_index[0] = MuonResiduals5DOFFitter::kResid;
1705  } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
1706  nvar = 1;
1707  vars_index[0] = MuonResiduals5DOFFitter::kResSlope;
1708  }
1709  } else if (fitter->type() == MuonResidualsFitter::k6DOF) {
1710  if (fitter->useRes() == MuonResidualsFitter::k1111) {
1711  nvar = 4;
1712  vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1713  vars_index[1] = MuonResiduals6DOFFitter::kResidY;
1714  vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
1715  vars_index[3] = MuonResiduals6DOFFitter::kResSlopeY;
1716  } else if (fitter->useRes() == MuonResidualsFitter::k1110) {
1717  nvar = 3;
1718  vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1719  vars_index[1] = MuonResiduals6DOFFitter::kResidY;
1720  vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
1721  } else if (fitter->useRes() == MuonResidualsFitter::k1010) {
1722  nvar = 2;
1723  vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1724  vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
1725  } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
1726  nvar = 2;
1727  vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1728  vars_index[1] = MuonResiduals6DOFFitter::kResidY;
1729  } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
1730  nvar = 1;
1731  vars_index[0] = MuonResiduals6DOFFitter::kResSlopeX;
1732  }
1733  } else if (fitter->type() == MuonResidualsFitter::k6DOFrphi) {
1734  if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 ||
1735  fitter->useRes() == MuonResidualsFitter::k1010) {
1736  nvar = 2;
1737  vars_index[0] = MuonResiduals6DOFrphiFitter::kResid;
1738  vars_index[1] = MuonResiduals6DOFrphiFitter::kResSlope;
1739  } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
1740  nvar = 1;
1741  vars_index[0] = MuonResiduals6DOFrphiFitter::kResid;
1742  } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
1743  nvar = 1;
1744  vars_index[0] = MuonResiduals6DOFrphiFitter::kResSlope;
1745  }
1746  } else
1747  assert(false);
1748 
1749  if (m_debug)
1750  std::cout << "selecting in " << chamberPrettyNameFromId(*index) << std::endl;
1751 
1752  fitter->selectPeakResiduals(m_peakNSigma, nvar, vars_index);
1753  }
1754 }
1755 
1756 std::string MuonAlignmentFromReference::chamberPrettyNameFromId(unsigned int idx) {
1757  DetId id(idx);
1758  char cname[40];
1759  if (id.subdetId() == MuonSubdetId::DT) {
1760  DTChamberId chamberId(id.rawId());
1761  sprintf(cname, "MB%+d/%d/%02d", chamberId.wheel(), chamberId.station(), chamberId.sector());
1762  } else if (id.subdetId() == MuonSubdetId::CSC) {
1763  CSCDetId chamberId(id.rawId());
1764  sprintf(cname,
1765  "ME%s%d/%d/%02d",
1766  (chamberId.endcap() == 1 ? "+" : "-"),
1767  chamberId.station(),
1768  chamberId.ring(),
1769  chamberId.chamber());
1770  }
1771  return std::string(cname);
1772 }
1773 
1774 void MuonAlignmentFromReference::fillNtuple() {
1775  // WARNING: does not support two bin option!!!
1776 
1777  for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1778  DetId detid(*index);
1779  if (detid.det() != DetId::Muon || !(detid.subdetId() == MuonSubdetId::DT || detid.subdetId() == MuonSubdetId::CSC))
1780  assert(false);
1781 
1782  if (detid.subdetId() == MuonSubdetId::DT) {
1783  m_tree_row.is_dt = (Bool_t) true;
1784  DTChamberId id(*index);
1785  m_tree_row.is_plus = (Bool_t) true;
1786  m_tree_row.station = (UChar_t)id.station();
1787  m_tree_row.ring_wheel = (Char_t)id.wheel();
1788  m_tree_row.sector = (UChar_t)id.sector();
1789  } else {
1790  m_tree_row.is_dt = (Bool_t) false;
1791  CSCDetId id(*index);
1792  m_tree_row.is_plus = (Bool_t)(id.endcap() == 1);
1793  m_tree_row.station = (UChar_t)id.station();
1794  m_tree_row.ring_wheel = (Char_t)id.ring();
1795  m_tree_row.sector = (UChar_t)id.chamber();
1796  }
1797 
1798  MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1799 
1800  std::vector<double*>::const_iterator residual = fitter->residualsPos_begin();
1801  std::vector<bool>::const_iterator residual_ok = fitter->residualsPos_ok_begin();
1802  for (; residual != fitter->residualsPos_end(); ++residual, ++residual_ok) {
1803  if (fitter->type() == MuonResidualsFitter::k5DOF || fitter->type() == MuonResidualsFitter::k6DOFrphi) {
1804  m_tree_row.res_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kResid];
1805  m_tree_row.res_y = (Float_t)0.;
1806  m_tree_row.res_slope_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kResSlope];
1807  m_tree_row.res_slope_y = (Float_t)0.;
1808  m_tree_row.pos_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPositionX];
1809  m_tree_row.pos_y = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPositionY];
1810  m_tree_row.angle_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kAngleX];
1811  m_tree_row.angle_y = (Float_t)(*residual)[MuonResiduals5DOFFitter::kAngleY];
1812  m_tree_row.pz = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPz];
1813  m_tree_row.pt = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPt];
1814  m_tree_row.q = (Char_t)(*residual)[MuonResiduals5DOFFitter::kCharge];
1815  m_tree_row.select = (Bool_t)*residual_ok;
1816  } else if (fitter->type() == MuonResidualsFitter::k6DOF) {
1817  m_tree_row.res_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResidX];
1818  m_tree_row.res_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResidY];
1819  m_tree_row.res_slope_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResSlopeX];
1820  m_tree_row.res_slope_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResSlopeY];
1821  m_tree_row.pos_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPositionX];
1822  m_tree_row.pos_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPositionY];
1823  m_tree_row.angle_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kAngleX];
1824  m_tree_row.angle_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kAngleY];
1825  m_tree_row.pz = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPz];
1826  m_tree_row.pt = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPt];
1827  m_tree_row.q = (Char_t)(*residual)[MuonResiduals6DOFFitter::kCharge];
1828  m_tree_row.select = (Bool_t)*residual_ok;
1829  } else
1830  assert(false);
1831 
1832  m_ttree->Fill();
1833  }
1834  }
1835 }
1836 
1837 void MuonAlignmentFromReference::parseReference(align::Alignables& reference,
1838  const align::Alignables& all_DT_chambers,
1839  const align::Alignables& all_CSC_chambers) {
1840  std::map<Alignable*, bool> already_seen;
1841 
1842  for (std::vector<std::string>::const_iterator name = m_reference.begin(); name != m_reference.end(); ++name) {
1843  bool parsing_error = false;
1844 
1845  bool barrel = (name->substr(0, 2) == std::string("MB"));
1846  bool endcap = (name->substr(0, 2) == std::string("ME"));
1847  if (!barrel && !endcap)
1848  parsing_error = true;
1849 
1850  if (!parsing_error && barrel) {
1851  int index = 2;
1852  if (name->substr(index, 1) == std::string(" "))
1853  index++;
1854 
1855  bool plus = true;
1856  if (name->substr(index, 1) == std::string("+")) {
1857  plus = true;
1858  index++;
1859  } else if (name->substr(index, 1) == std::string("-")) {
1860  plus = false;
1861  index++;
1862  } else if (numeric(name->substr(index, 1))) {
1863  } else
1864  parsing_error = true;
1865 
1866  int wheel = 0;
1867  bool wheel_digit = false;
1868  while (!parsing_error && numeric(name->substr(index, 1))) {
1869  wheel *= 10;
1870  wheel += number(name->substr(index, 1));
1871  wheel_digit = true;
1872  index++;
1873  }
1874  if (!plus)
1875  wheel *= -1;
1876  if (!wheel_digit)
1877  parsing_error = true;
1878 
1879  if (name->substr(index, 1) != std::string(" "))
1880  parsing_error = true;
1881  index++;
1882 
1883  int station = 0;
1884  bool station_digit = false;
1885  while (!parsing_error && numeric(name->substr(index, 1))) {
1886  station *= 10;
1887  station += number(name->substr(index, 1));
1888  station_digit = true;
1889  index++;
1890  }
1891  if (!station_digit)
1892  parsing_error = true;
1893 
1894  if (name->substr(index, 1) != std::string(" "))
1895  parsing_error = true;
1896  index++;
1897 
1898  int sector = 0;
1899  bool sector_digit = false;
1900  while (!parsing_error && numeric(name->substr(index, 1))) {
1901  sector *= 10;
1902  sector += number(name->substr(index, 1));
1903  sector_digit = true;
1904  index++;
1905  }
1906  if (!sector_digit)
1907  parsing_error = true;
1908 
1909  if (!parsing_error) {
1910  bool no_such_chamber = false;
1911 
1912  if (wheel < -2 || wheel > 2)
1913  no_such_chamber = true;
1914  if (station < 1 || station > 4)
1915  no_such_chamber = true;
1916  if (station == 4 && (sector < 1 || sector > 14))
1917  no_such_chamber = true;
1918  if (station < 4 && (sector < 1 || sector > 12))
1919  no_such_chamber = true;
1920 
1921  if (no_such_chamber)
1922  throw cms::Exception("MuonAlignmentFromReference")
1923  << "reference chamber doesn't exist: " << (*name) << std::endl;
1924 
1925  DTChamberId id(wheel, station, sector);
1926  for (const auto& ali : all_DT_chambers) {
1927  if (ali->geomDetId().rawId() == id.rawId()) {
1928  std::map<Alignable*, bool>::const_iterator trial = already_seen.find(ali);
1929  if (trial == already_seen.end()) {
1930  reference.push_back(ali);
1931  already_seen[ali] = true;
1932  }
1933  }
1934  }
1935  } // if (!parsing_error)
1936  }
1937 
1938  if (!parsing_error && endcap) {
1939  int index = 2;
1940  if (name->substr(index, 1) == std::string(" "))
1941  index++;
1942 
1943  bool plus = true;
1944  if (name->substr(index, 1) == std::string("+")) {
1945  plus = true;
1946  index++;
1947  } else if (name->substr(index, 1) == std::string("-")) {
1948  plus = false;
1949  index++;
1950  } else if (numeric(name->substr(index, 1))) {
1951  } else
1952  parsing_error = true;
1953 
1954  int station = 0;
1955  bool station_digit = false;
1956  while (!parsing_error && numeric(name->substr(index, 1))) {
1957  station *= 10;
1958  station += number(name->substr(index, 1));
1959  station_digit = true;
1960  index++;
1961  }
1962  if (!plus)
1963  station *= -1;
1964  if (!station_digit)
1965  parsing_error = true;
1966 
1967  if (name->substr(index, 1) != std::string("/"))
1968  parsing_error = true;
1969  index++;
1970 
1971  int ring = 0;
1972  bool ring_digit = false;
1973  while (!parsing_error && numeric(name->substr(index, 1))) {
1974  ring *= 10;
1975  ring += number(name->substr(index, 1));
1976  ring_digit = true;
1977  index++;
1978  }
1979  if (!ring_digit)
1980  parsing_error = true;
1981 
1982  if (name->substr(index, 1) != std::string(" "))
1983  parsing_error = true;
1984  index++;
1985 
1986  int chamber = 0;
1987  bool chamber_digit = false;
1988  while (!parsing_error && numeric(name->substr(index, 1))) {
1989  chamber *= 10;
1990  chamber += number(name->substr(index, 1));
1991  chamber_digit = true;
1992  index++;
1993  }
1994  if (!chamber_digit)
1995  parsing_error = true;
1996 
1997  if (!parsing_error) {
1998  bool no_such_chamber = false;
1999 
2000  int endcap = (station > 0 ? 1 : 2);
2001  station = abs(station);
2002  if (station < 1 || station > 4)
2003  no_such_chamber = true;
2004  if (station == 1 && (ring < 1 || ring > 4))
2005  no_such_chamber = true;
2006  if (station > 1 && (ring < 1 || ring > 2))
2007  no_such_chamber = true;
2008  if (station == 1 && (chamber < 1 || chamber > 36))
2009  no_such_chamber = true;
2010  if (station > 1 && ring == 1 && (chamber < 1 || chamber > 18))
2011  no_such_chamber = true;
2012  if (station > 1 && ring == 2 && (chamber < 1 || chamber > 36))
2013  no_such_chamber = true;
2014 
2015  if (no_such_chamber)
2016  throw cms::Exception("MuonAlignmentFromReference")
2017  << "reference chamber doesn't exist: " << (*name) << std::endl;
2018 
2019  CSCDetId id(endcap, station, ring, chamber);
2020  for (const auto& ali : all_CSC_chambers) {
2021  if (ali->geomDetId().rawId() == id.rawId()) {
2022  std::map<Alignable*, bool>::const_iterator trial = already_seen.find(ali);
2023  if (trial == already_seen.end()) {
2024  reference.push_back(ali);
2025  already_seen[ali] = true;
2026  }
2027  }
2028  }
2029  } // if (!parsing_error)
2030  } // endcap
2031 
2032  if (parsing_error)
2033  throw cms::Exception("MuonAlignmentFromReference")
2034  << "reference chamber name is malformed: " << (*name) << std::endl;
2035  }
2036 }
2037 
2038 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
2039 DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory, MuonAlignmentFromReference, "MuonAlignmentFromReference");
MuonAlignmentFromReference::initialize
void initialize(const edm::EventSetup &iSetup, AlignableTracker *alignableTracker, AlignableMuon *alignableMuon, AlignableExtras *extras, AlignmentParameterStore *alignmentParameterStore) override
Call at beginning of job (must be implemented in derived class)
Definition: MuonAlignmentFromReference.cc:308
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
CSCDetId::chamber
int chamber() const
Definition: CSCDetId.h:62
MuonResidualsTwoBin::residualsPos_end
std::vector< double * >::const_iterator residualsPos_end() const
Definition: MuonResidualsTwoBin.h:257
DTChamberId
DTChamberId
Definition: RecoLocalMuon_RECO.doi:23
MuonResidualsTwoBin::residualsPos_begin
std::vector< double * >::const_iterator residualsPos_begin() const
Definition: MuonResidualsTwoBin.h:256
looper.cfg
tuple cfg
Definition: looper.py:296
mps_fire.i
i
Definition: mps_fire.py:428
MuonAlignmentFromReference::m_MagFieldToken
const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > m_MagFieldToken
Definition: MuonAlignmentFromReference.cc:116
DiDispStaMuonMonitor_cfi.pt
tuple pt
Definition: DiDispStaMuonMonitor_cfi.py:39
gpuClustering::id
uint16_t *__restrict__ id
Definition: gpuClusterChargeCut.h:20
Reference_intrackfit_cff.barrel
barrel
Definition: Reference_intrackfit_cff.py:37
DetIdAssociatorESProducer_cff.muonDetIdAssociator
tuple muonDetIdAssociator
Definition: DetIdAssociatorESProducer_cff.py:34
MuonAlignmentFromReference::m_DetIdToken
const edm::ESGetToken< DetIdAssociator, DetIdAssociatorRecord > m_DetIdToken
Definition: MuonAlignmentFromReference.cc:118
AlignableMuon::CSCChambers
align::Alignables CSCChambers()
Definition: AlignableMuon.cc:418
Reference_intrackfit_cff.endcap
endcap
Definition: Reference_intrackfit_cff.py:82
AlignmentParameterStore::setAlignmentPositionError
void setAlignmentPositionError(const align::Alignables &alivec, double valshift, double valrot)
Set Alignment position error.
Definition: AlignmentParameterStore.cc:646
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
HLT_FULL_cff.magneticField
tuple magneticField
Definition: HLT_FULL_cff.py:360
MuonChamberResidual::ChambW
double ChambW() const
Definition: MuonChamberResidual.h:67
MuonResidualsTwoBin::residualsPos_ok_begin
std::vector< bool >::const_iterator residualsPos_ok_begin() const
Definition: MuonResidualsTwoBin.h:261
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
AlignmentParameterStore::applyParameters
void applyParameters(void)
Obsolete: Use AlignableNavigator::alignableDetFromDetId and alignableFromAlignableDet.
Definition: AlignmentParameterStore.cc:255
HLT_FULL_cff.track
tuple track
Definition: HLT_FULL_cff.py:11953
cms::cuda::assert
assert(be >=bs)
AlignmentAlgorithmBase.h
Interface/Base class for alignment algorithms, each alignment algorithm has to be derived from this c...
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
MuonChamberResidual::trackx
double trackx() const
Definition: MuonChamberResidual.h:64
AlignmentAlgorithmBase::EventInfo::trajTrackPairs
const ConstTrajTrackPairCollection & trajTrackPairs() const
Definition: AlignmentAlgorithmBase.h:85
AlignmentAlgorithmBase::EventInfo
define event information passed to algorithms
Definition: AlignmentAlgorithmBase.h:73
MuonAlignmentFromReference::m_alignmentParameterStore
AlignmentParameterStore * m_alignmentParameterStore
Definition: MuonAlignmentFromReference.cc:156
MuonResidualsFromTrack::chamberIds
const std::vector< DetId > chamberIds() const
Definition: MuonResidualsFromTrack.h:89
CSCDetId::endcap
int endcap() const
Definition: CSCDetId.h:85
MuonAlignmentFromReference::processMuonResidualsFromTrack
void processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft)
Definition: MuonAlignmentFromReference.cc:529
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:128
zeeHLT_cff.report
tuple report
Definition: zeeHLT_cff.py:9
MuonAlignmentFromReference::m_tree_row
MuonResidualsFitter::MuonAlignmentTreeRow m_tree_row
Definition: MuonAlignmentFromReference.cc:191
MuonChamberResidual::resslope
double resslope() const
Definition: MuonChamberResidual.h:56
align_cfg.residualsModel
list residualsModel
Definition: align_cfg.py:34
MuonChamberResidual::trackdxdz
double trackdxdz() const
Definition: MuonChamberResidual.h:62
MuonAlignmentFromReference::m_indexes
std::vector< unsigned int > m_indexes
Definition: MuonAlignmentFromReference.cc:160
RecoTauValidation_cfi.reference
tuple reference
Definition: RecoTauValidation_cfi.py:233
MuonAlignmentFromReference::chamberPrettyNameFromId
std::string chamberPrettyNameFromId(unsigned int idx)
Definition: MuonAlignmentFromReference.cc:1756
AlignableMuon::DTChambers
align::Alignables DTChambers()
Definition: AlignableMuon.cc:375
MillePedeFileConverter_cfg.fileName
tuple fileName
Definition: MillePedeFileConverter_cfg.py:32
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:637
MuonResidualsFromTrack::getTrack
const reco::Track * getTrack()
Definition: MuonResidualsFromTrack.h:78
MuonAlignmentFromReference::m_builderToken
const MuonResidualsFromTrack::BuilderToken m_builderToken
Definition: MuonAlignmentFromReference.cc:119
MuonChamberResidual::chamberAlignable
AlignableDetOrUnitPtr chamberAlignable() const
Definition: MuonChamberResidual.h:50
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::Transition
Transition
Definition: Transition.h:12
MuonAlignmentFromReference::m_propToken
const edm::ESGetToken< Propagator, TrackingComponentsRecord > m_propToken
Definition: MuonAlignmentFromReference.cc:117
MuonAlignmentFromReference::parseReference
void parseReference(align::Alignables &reference, const align::Alignables &all_DT_chambers, const align::Alignables &all_CSC_chambers)
Definition: MuonAlignmentFromReference.cc:1837
MuonResidualsTwoBin::write
void write(FILE *file, int which=0)
Definition: MuonResidualsTwoBin.h:102
MuonAlignmentFromReference::MuonAlignmentFromReference
MuonAlignmentFromReference(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC)
Definition: MuonAlignmentFromReference.cc:196
MuonAlignmentFromReference::m_cscGeometryToken
const edm::ESGetToken< CSCGeometry, MuonGeometryRecord > m_cscGeometryToken
Definition: MuonAlignmentFromReference.cc:114
MuonChamberResidual::tracky
double tracky() const
Definition: MuonChamberResidual.h:65
MuonAlignmentFromReference::m_fitterOrder
std::map< unsigned int, MuonResidualsTwoBin * > m_fitterOrder
Definition: MuonAlignmentFromReference.cc:161
MuonAlignmentFromReference::m_reference
std::vector< std::string > m_reference
Definition: MuonAlignmentFromReference.cc:123
MuonAlignmentFromReference::m_fitters
std::map< Alignable *, MuonResidualsTwoBin * > m_fitters
Definition: MuonAlignmentFromReference.cc:159
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:646
CSCDetId::ring
int ring() const
Definition: CSCDetId.h:68
TFileService::file
TFile & file() const
return opened TFile
Definition: TFileService.h:37
DetId
Definition: DetId.h:17
MuonAlignmentFromReference::run
void run(const edm::EventSetup &iSetup, const EventInfo &eventInfo) override
Run the algorithm (must be implemented in derived class)
Definition: MuonAlignmentFromReference.cc:441
AlgebraicVector
CLHEP::HepVector AlgebraicVector
Definition: AlgebraicObjects.h:13
TFileService::mkdir
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
align_cfg.useResiduals
list useResiduals
Definition: align_cfg.py:36
MuonChamberResidual::trackdydz
double trackdydz() const
Definition: MuonChamberResidual.h:63
MuonResidualsTwoBin::read
void read(FILE *file, int which=0)
Definition: MuonResidualsTwoBin.h:110
align::Alignables
std::vector< Alignable * > Alignables
Definition: Utilities.h:31
MuonAlignmentFromReference::m_me11map
std::map< Alignable *, Alignable * > m_me11map
Definition: MuonAlignmentFromReference.cc:158
edm::InputTag::label
std::string const & label() const
Definition: InputTag.h:36
HLT_FULL_cff.InputTag
tuple InputTag
Definition: HLT_FULL_cff.py:86956
MuonAlignmentFromReference::m_globTackingToken
const edm::ESGetToken< GlobalTrackingGeometry, GlobalTrackingGeometryRecord > m_globTackingToken
Definition: MuonAlignmentFromReference.cc:115
AlignmentAlgorithmBase::terminate
virtual void terminate()
Called at end of job (must be implemented in derived class)
Definition: AlignmentAlgorithmBase.h:158
DTChamberId::sector
int sector() const
Definition: DTChamberId.h:49
AlignmentAlgorithmBase::EventInfo::beamSpot
const reco::BeamSpot & beamSpot() const
Definition: AlignmentAlgorithmBase.h:86
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:15
MuonAlignmentFromReference::m_readTemporaryFiles
std::vector< std::string > m_readTemporaryFiles
Definition: MuonAlignmentFromReference.cc:137
CSCDetId::station
int station() const
Definition: CSCDetId.h:79
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
MuonChamberResidual::Chambl
double Chambl() const
Definition: MuonChamberResidual.h:68
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:596
DEFINE_EDM_PLUGIN
#define DEFINE_EDM_PLUGIN(factory, type, name)
Definition: PluginFactory.h:124
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:59
compare.tfile
tuple tfile
Definition: compare.py:324
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:157
MuonSubdetId::DT
static constexpr int DT
Definition: MuonSubdetId.h:11
DTChamberId::station
int station() const
Return the station number.
Definition: DTChamberId.h:42
MuonResidualsTwoBin::selectPeakResiduals
void selectPeakResiduals(double nsigma, int nvar, int *vars)
Definition: MuonResidualsTwoBin.h:231
AlignmentParameters::cloneFromSelected
virtual AlignmentParameters * cloneFromSelected(const AlgebraicVector &par, const AlgebraicSymMatrix &cov) const =0
DTChamberId::wheel
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:39
reco::TrackBase::dxy
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:608
AlignableMuon
Constructor of the full muon geometry.
Definition: AlignableMuon.h:38
MuonSubdetId::CSC
static constexpr int CSC
Definition: MuonSubdetId.h:12
MuonChamberResidual::residual
double residual() const
Definition: MuonChamberResidual.h:54
DeDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
Alignable::geomDetId
const DetId & geomDetId() const
Definition: Alignable.h:177
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:443
AlignmentParameterStore::alignables
const align::Alignables & alignables(void) const
get all alignables
Definition: AlignmentParameterStore.h:47
DetId::det
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
MuonResidualsFromTrack::chamberResidual
MuonChamberResidual * chamberResidual(DetId chamberId, int type)
Definition: MuonResidualsFromTrack.cc:735
AlignmentAlgorithmBase::ConstTrajTrackPairCollection
std::vector< ConstTrajTrackPair > ConstTrajTrackPairCollection
Definition: AlignmentAlgorithmBase.h:68
Generated for CMSSW Reference Manual by doxygen 1.8.5