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