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