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