Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:34:34 

0001 // -*- C++ -*-
0002 //
0003 // Package:    MuonAlignmentAlgorithms
0004 // Class:      MuonAlignmentFromReference
0005 //
0006 /**\class MuonAlignmentFromReference MuonAlignmentFromReference.cc Alignment/MuonAlignmentFromReference/interface/MuonAlignmentFromReference.h
0007 
0008 Description: <one line class summary>
0009 
0010 Implementation:
0011 <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Jim Pivarski,,,
0015 //         Created:  Sat Jan 24 16:20:28 CST 2009
0016 // $Id: MuonAlignmentFromReference.cc,v 1.39 2011/10/13 00:03:12 khotilov Exp $
0017 
0018 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmBase.h"
0019 
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0027 
0028 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0029 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0030 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0031 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0032 
0033 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0034 #include "MagneticField/Engine/interface/MagneticField.h"
0035 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0036 #include "TrackingTools/Records/interface/DetIdAssociatorRecord.h"
0037 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0038 
0039 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0040 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0041 #include "DataFormats/MuonDetId/interface/DTSuperLayerId.h"
0042 #include "DataFormats/TrackReco/interface/Track.h"
0043 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0044 
0045 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0046 
0047 #include "Alignment/CommonAlignment/interface/Alignable.h"
0048 #include "Alignment/CommonAlignment/interface/AlignableDetUnit.h"
0049 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
0050 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
0051 #include "Alignment/MuonAlignment/interface/AlignableDTSuperLayer.h"
0052 #include "Alignment/MuonAlignment/interface/AlignableDTChamber.h"
0053 #include "Alignment/MuonAlignment/interface/AlignableDTStation.h"
0054 #include "Alignment/MuonAlignment/interface/AlignableDTWheel.h"
0055 #include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
0056 #include "Alignment/MuonAlignment/interface/AlignableCSCRing.h"
0057 #include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"
0058 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
0059 
0060 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
0061 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
0062 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals5DOFFitter.h"
0063 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFrphiFitter.h"
0064 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h"
0065 
0066 #include "TFile.h"
0067 #include "TTree.h"
0068 #include "TStopwatch.h"
0069 
0070 #include <map>
0071 #include <sstream>
0072 #include <fstream>
0073 
0074 class MuonAlignmentFromReference : public AlignmentAlgorithmBase {
0075 public:
0076   MuonAlignmentFromReference(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC);
0077   ~MuonAlignmentFromReference() override;
0078 
0079   void initialize(const edm::EventSetup& iSetup,
0080                   AlignableTracker* alignableTracker,
0081                   AlignableMuon* alignableMuon,
0082                   AlignableExtras* extras,
0083                   AlignmentParameterStore* alignmentParameterStore) override;
0084 
0085   void startNewLoop() override {}
0086 
0087   void run(const edm::EventSetup& iSetup, const EventInfo& eventInfo) override;
0088 
0089   void processMuonResidualsFromTrack(MuonResidualsFromTrack& mrft);
0090 
0091   void terminate(const edm::EventSetup& iSetup) override;
0092 
0093 private:
0094   bool numeric(std::string s);
0095   int number(std::string s);
0096   std::string chamberPrettyNameFromId(unsigned int idx);
0097 
0098   void parseReference(align::Alignables& reference,
0099                       const align::Alignables& all_DT_chambers,
0100                       const align::Alignables& all_CSC_chambers);
0101 
0102   void fitAndAlign();
0103   void readTmpFiles();
0104   void writeTmpFiles();
0105 
0106   void selectResidualsPeaks();
0107   void correctBField();
0108   void fiducialCuts();
0109   void eraseNotSelectedResiduals();
0110 
0111   void fillNtuple();
0112 
0113   // tokens
0114   const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> m_cscGeometryToken;
0115   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_globTackingToken;
0116   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagFieldToken;
0117   const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_propToken;
0118   const edm::ESGetToken<DetIdAssociator, DetIdAssociatorRecord> m_DetIdToken;
0119   const MuonResidualsFromTrack::BuilderToken m_builderToken;
0120 
0121   // configutarion paramenters:
0122   edm::InputTag m_muonCollectionTag;
0123   std::vector<std::string> m_reference;
0124   double m_minTrackPt;
0125   double m_maxTrackPt;
0126   double m_minTrackP;
0127   double m_maxTrackP;
0128   double m_maxDxy;
0129   int m_minTrackerHits;
0130   double m_maxTrackerRedChi2;
0131   bool m_allowTIDTEC;
0132   int m_minNCrossedChambers;
0133   int m_minDT13Hits;
0134   int m_minDT2Hits;
0135   int m_minCSCHits;
0136   std::string m_writeTemporaryFile;
0137   std::vector<std::string> m_readTemporaryFiles;
0138   bool m_doAlignment;
0139   int m_strategy;
0140   std::string m_residualsModel;
0141   int m_minAlignmentHits;
0142   bool m_twoBin;
0143   bool m_combineME11;
0144   bool m_weightAlignment;
0145   std::string m_reportFileName;
0146   double m_maxResSlopeY;
0147   bool m_createNtuple;
0148   double m_peakNSigma;
0149   int m_BFieldCorrection;
0150   bool m_doDT;
0151   bool m_doCSC;
0152   std::string m_useResiduals;
0153 
0154   // utility objects
0155   AlignableNavigator* m_alignableNavigator;
0156   AlignmentParameterStore* m_alignmentParameterStore;
0157   align::Alignables m_alignables;
0158   std::map<Alignable*, Alignable*> m_me11map;
0159   std::map<Alignable*, MuonResidualsTwoBin*> m_fitters;
0160   std::vector<unsigned int> m_indexes;
0161   std::map<unsigned int, MuonResidualsTwoBin*> m_fitterOrder;
0162 
0163   // counters
0164   long m_counter_events;
0165   long m_counter_tracks;
0166   long m_counter_trackmomentum;
0167   long m_counter_trackdxy;
0168   long m_counter_trackerhits;
0169   long m_counter_trackerchi2;
0170   long m_counter_trackertidtec;
0171   long m_counter_minchambers;
0172   long m_counter_totchambers;
0173   long m_counter_station123;
0174   long m_counter_station123valid;
0175   long m_counter_station123dt13hits;
0176   long m_counter_station123dt2hits;
0177   long m_counter_station123aligning;
0178   long m_counter_station4;
0179   long m_counter_station4valid;
0180   long m_counter_station4hits;
0181   long m_counter_station4aligning;
0182   long m_counter_csc;
0183   long m_counter_cscvalid;
0184   long m_counter_cschits;
0185   long m_counter_cscaligning;
0186   long m_counter_resslopey;
0187 
0188   // debug ntuple
0189   void bookNtuple();
0190   TTree* m_ttree;
0191   MuonResidualsFitter::MuonAlignmentTreeRow m_tree_row;
0192 
0193   bool m_debug;
0194 };
0195 
0196 MuonAlignmentFromReference::MuonAlignmentFromReference(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0197     : AlignmentAlgorithmBase(cfg, iC),
0198       m_cscGeometryToken(iC.esConsumes<edm::Transition::BeginRun>()),
0199       m_globTackingToken(iC.esConsumes()),
0200       m_MagFieldToken(iC.esConsumes()),
0201       m_propToken(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0202       m_DetIdToken(iC.esConsumes(edm::ESInputTag("", "MuonDetIdAssociator"))),
0203       m_builderToken(iC.esConsumes(MuonResidualsFromTrack::builderESInputTag())),
0204       m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag")),
0205       m_reference(cfg.getParameter<std::vector<std::string> >("reference")),
0206       m_minTrackPt(cfg.getParameter<double>("minTrackPt")),
0207       m_maxTrackPt(cfg.getParameter<double>("maxTrackPt")),
0208       m_minTrackP(cfg.getParameter<double>("minTrackP")),
0209       m_maxTrackP(cfg.getParameter<double>("maxTrackP")),
0210       m_maxDxy(cfg.getParameter<double>("maxDxy")),
0211       m_minTrackerHits(cfg.getParameter<int>("minTrackerHits")),
0212       m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2")),
0213       m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC")),
0214       m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers")),
0215       m_minDT13Hits(cfg.getParameter<int>("minDT13Hits")),
0216       m_minDT2Hits(cfg.getParameter<int>("minDT2Hits")),
0217       m_minCSCHits(cfg.getParameter<int>("minCSCHits")),
0218       m_writeTemporaryFile(cfg.getParameter<std::string>("writeTemporaryFile")),
0219       m_readTemporaryFiles(cfg.getParameter<std::vector<std::string> >("readTemporaryFiles")),
0220       m_doAlignment(cfg.getParameter<bool>("doAlignment")),
0221       m_strategy(cfg.getParameter<int>("strategy")),
0222       m_residualsModel(cfg.getParameter<std::string>("residualsModel")),
0223       m_minAlignmentHits(cfg.getParameter<int>("minAlignmentHits")),
0224       m_twoBin(cfg.getParameter<bool>("twoBin")),
0225       m_combineME11(cfg.getParameter<bool>("combineME11")),
0226       m_weightAlignment(cfg.getParameter<bool>("weightAlignment")),
0227       m_reportFileName(cfg.getParameter<std::string>("reportFileName")),
0228       m_maxResSlopeY(cfg.getParameter<double>("maxResSlopeY")),
0229       m_createNtuple(cfg.getParameter<bool>("createNtuple")),
0230       m_peakNSigma(cfg.getParameter<double>("peakNSigma")),
0231       m_BFieldCorrection(cfg.getParameter<int>("bFieldCorrection")),
0232       m_doDT(cfg.getParameter<bool>("doDT")),
0233       m_doCSC(cfg.getParameter<bool>("doCSC")),
0234       m_useResiduals(cfg.getParameter<std::string>("useResiduals")) {
0235   // alignment requires a TFile to provide plots to check the fit output
0236   // just filling the residuals lists does not
0237   // but we don't want to wait until the end of the job to find out that the TFile is missing
0238   if (m_doAlignment || m_createNtuple) {
0239     edm::Service<TFileService> fs;
0240     TFile& tfile = fs->file();
0241     tfile.ls();
0242   }
0243 
0244   m_ttree = nullptr;
0245   if (m_createNtuple)
0246     bookNtuple();
0247 
0248   m_counter_events = 0;
0249   m_counter_tracks = 0;
0250   m_counter_trackmomentum = 0;
0251   m_counter_trackdxy = 0;
0252   m_counter_trackerhits = 0;
0253   m_counter_trackerchi2 = 0;
0254   m_counter_trackertidtec = 0;
0255   m_counter_minchambers = 0;
0256   m_counter_totchambers = 0;
0257   m_counter_station123 = 0;
0258   m_counter_station123valid = 0;
0259   m_counter_station123dt13hits = 0;
0260   m_counter_station123dt2hits = 0;
0261   m_counter_station123aligning = 0;
0262   m_counter_station4 = 0;
0263   m_counter_station4valid = 0;
0264   m_counter_station4hits = 0;
0265   m_counter_station4aligning = 0;
0266   m_counter_csc = 0;
0267   m_counter_cscvalid = 0;
0268   m_counter_cschits = 0;
0269   m_counter_cscaligning = 0;
0270   m_counter_resslopey = 0;
0271 
0272   m_debug = false;
0273 }
0274 
0275 MuonAlignmentFromReference::~MuonAlignmentFromReference() { delete m_alignableNavigator; }
0276 
0277 void MuonAlignmentFromReference::bookNtuple() {
0278   edm::Service<TFileService> fs;
0279   m_ttree = fs->make<TTree>("mual_ttree", "mual_ttree");
0280   m_ttree->Branch("is_plus", &m_tree_row.is_plus, "is_plus/O");
0281   m_ttree->Branch("is_dt", &m_tree_row.is_dt, "is_dt/O");
0282   m_ttree->Branch("station", &m_tree_row.station, "station/b");
0283   m_ttree->Branch("ring_wheel", &m_tree_row.ring_wheel, "ring_wheel/B");
0284   m_ttree->Branch("sector", &m_tree_row.sector, "sector/b");
0285   m_ttree->Branch("res_x", &m_tree_row.res_x, "res_x/F");
0286   m_ttree->Branch("res_y", &m_tree_row.res_y, "res_y/F");
0287   m_ttree->Branch("res_slope_x", &m_tree_row.res_slope_x, "res_slope_x/F");
0288   m_ttree->Branch("res_slope_y", &m_tree_row.res_slope_y, "res_slope_y/F");
0289   m_ttree->Branch("pos_x", &m_tree_row.pos_x, "pos_x/F");
0290   m_ttree->Branch("pos_y", &m_tree_row.pos_y, "pos_y/F");
0291   m_ttree->Branch("angle_x", &m_tree_row.angle_x, "angle_x/F");
0292   m_ttree->Branch("angle_y", &m_tree_row.angle_y, "angle_y/F");
0293   m_ttree->Branch("pz", &m_tree_row.pz, "pz/F");
0294   m_ttree->Branch("pt", &m_tree_row.pt, "pt/F");
0295   m_ttree->Branch("q", &m_tree_row.q, "q/B");
0296   m_ttree->Branch("select", &m_tree_row.select, "select/O");
0297   //m_ttree->Branch("",&m_tree_row.,"/");
0298 }
0299 
0300 bool MuonAlignmentFromReference::numeric(std::string s) { return s.length() == 1 && std::isdigit(s[0]); }
0301 
0302 int MuonAlignmentFromReference::number(std::string s) {
0303   if (!numeric(s))
0304     assert(false);
0305   return atoi(s.c_str());
0306 }
0307 
0308 void MuonAlignmentFromReference::initialize(const edm::EventSetup& iSetup,
0309                                             AlignableTracker* alignableTracker,
0310                                             AlignableMuon* alignableMuon,
0311                                             AlignableExtras* extras,
0312                                             AlignmentParameterStore* alignmentParameterStore) {
0313   if (alignableMuon == nullptr)
0314     throw cms::Exception("MuonAlignmentFromReference") << "doMuon must be set to True" << std::endl;
0315 
0316   m_alignableNavigator = new AlignableNavigator(alignableMuon);
0317   m_alignmentParameterStore = alignmentParameterStore;
0318   m_alignables = m_alignmentParameterStore->alignables();
0319 
0320   int residualsModel;
0321   if (m_residualsModel == std::string("pureGaussian"))
0322     residualsModel = MuonResidualsFitter::kPureGaussian;
0323   else if (m_residualsModel == std::string("pureGaussian2D"))
0324     residualsModel = MuonResidualsFitter::kPureGaussian2D;
0325   else if (m_residualsModel == std::string("powerLawTails"))
0326     residualsModel = MuonResidualsFitter::kPowerLawTails;
0327   else if (m_residualsModel == std::string("ROOTVoigt"))
0328     residualsModel = MuonResidualsFitter::kROOTVoigt;
0329   else if (m_residualsModel == std::string("GaussPowerTails"))
0330     residualsModel = MuonResidualsFitter::kGaussPowerTails;
0331   else
0332     throw cms::Exception("MuonAlignmentFromReference")
0333         << "unrecognized residualsModel: \"" << m_residualsModel << "\"" << std::endl;
0334 
0335   int useResiduals;
0336   if (m_useResiduals == std::string("1111"))
0337     useResiduals = MuonResidualsFitter::k1111;
0338   else if (m_useResiduals == std::string("1110"))
0339     useResiduals = MuonResidualsFitter::k1110;
0340   else if (m_useResiduals == std::string("1100"))
0341     useResiduals = MuonResidualsFitter::k1100;
0342   else if (m_useResiduals == std::string("1000"))
0343     useResiduals = MuonResidualsFitter::k1000;
0344   else if (m_useResiduals == std::string("1010"))
0345     useResiduals = MuonResidualsFitter::k1010;
0346   else if (m_useResiduals == std::string("0010"))
0347     useResiduals = MuonResidualsFitter::k0010;
0348   else
0349     throw cms::Exception("MuonAlignmentFromReference")
0350         << "unrecognized useResiduals: \"" << m_useResiduals << "\"" << std::endl;
0351 
0352   const CSCGeometry* cscGeometry = &iSetup.getData(m_cscGeometryToken);
0353 
0354   // set up the MuonResidualsFitters (which also collect residuals for fitting)
0355   m_me11map.clear();
0356   m_fitters.clear();
0357   m_indexes.clear();
0358   m_fitterOrder.clear();
0359 
0360   for (const auto& ali : m_alignables) {
0361     bool made_fitter = false;
0362 
0363     // fitters for DT
0364     if (ali->alignableObjectId() == align::AlignableDTChamber) {
0365       DTChamberId id(ali->geomDetId().rawId());
0366 
0367       if (id.station() == 4) {
0368         m_fitters[ali] = new MuonResidualsTwoBin(
0369             m_twoBin,
0370             new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment),
0371             new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment));
0372         made_fitter = true;
0373       } else {
0374         m_fitters[ali] = new MuonResidualsTwoBin(
0375             m_twoBin,
0376             new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment),
0377             new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment));
0378         made_fitter = true;
0379       }
0380     }
0381 
0382     // fitters for CSC
0383     else if (ali->alignableObjectId() == align::AlignableCSCChamber) {
0384       auto thisali = ali;
0385       CSCDetId id(ali->geomDetId().rawId());
0386 
0387       // take care of ME1/1a
0388       if (m_combineME11 && id.station() == 1 && id.ring() == 4) {
0389         CSCDetId pairid(id.endcap(), 1, 1, id.chamber());
0390 
0391         for (const auto& ali2 : m_alignables) {
0392           if (ali2->alignableObjectId() == align::AlignableCSCChamber && ali2->geomDetId().rawId() == pairid.rawId()) {
0393             thisali = ali2;
0394             break;
0395           }
0396         }
0397         m_me11map[ali] = thisali;  // points from each ME1/4 chamber to the corresponding ME1/1 chamber
0398       }
0399 
0400       if (thisali == ali)  // don't make fitters for ME1/4; they get taken care of in ME1/1
0401       {
0402         m_fitters[ali] = new MuonResidualsTwoBin(
0403             m_twoBin,
0404             new MuonResiduals6DOFrphiFitter(
0405                 residualsModel, m_minAlignmentHits, useResiduals, cscGeometry, m_weightAlignment),
0406             new MuonResiduals6DOFrphiFitter(
0407                 residualsModel, m_minAlignmentHits, useResiduals, cscGeometry, m_weightAlignment));
0408         made_fitter = true;
0409       }
0410     }
0411 
0412     else {
0413       throw cms::Exception("MuonAlignmentFromReference")
0414           << "only DTChambers and CSCChambers can be aligned with this module" << std::endl;
0415     }
0416 
0417     if (made_fitter) {
0418       m_fitters[ali]->setStrategy(m_strategy);
0419 
0420       int index = ali->geomDetId().rawId();
0421       m_indexes.push_back(index);
0422       m_fitterOrder[index] = m_fitters[ali];
0423     }
0424   }  // end loop over chambers chosen for alignment
0425 
0426   // cannonical order of fitters in the file
0427   std::sort(m_indexes.begin(), m_indexes.end());
0428 
0429   // de-weight all chambers but the reference
0430   const auto& all_DT_chambers = alignableMuon->DTChambers();
0431   const auto& all_CSC_chambers = alignableMuon->CSCChambers();
0432   align::Alignables reference;
0433   if (!m_reference.empty())
0434     parseReference(reference, all_DT_chambers, all_CSC_chambers);
0435 
0436   alignmentParameterStore->setAlignmentPositionError(all_DT_chambers, 100000000., 0.);
0437   alignmentParameterStore->setAlignmentPositionError(all_CSC_chambers, 100000000., 0.);
0438   alignmentParameterStore->setAlignmentPositionError(reference, 0., 0.);
0439 }
0440 
0441 void MuonAlignmentFromReference::run(const edm::EventSetup& iSetup, const EventInfo& eventInfo) {
0442   if (m_debug)
0443     std::cout << "****** EVENT START *******" << std::endl;
0444   m_counter_events++;
0445 
0446   const GlobalTrackingGeometry* globalGeometry = &iSetup.getData(m_globTackingToken);
0447   const MagneticField* magneticField = &iSetup.getData(m_MagFieldToken);
0448   const Propagator* prop = &iSetup.getData(m_propToken);
0449   const DetIdAssociator* muonDetIdAssociator = &iSetup.getData(m_DetIdToken);
0450   auto builder = iSetup.getHandle(m_builderToken);
0451 
0452   if (m_muonCollectionTag.label().empty())  // use trajectories
0453   {
0454     if (m_debug)
0455       std::cout << "JUST BEFORE LOOP OVER trajTrackPairs" << std::endl;
0456     // const ConstTrajTrackPairCollection &trajtracks = eventInfo.trajTrackPairs_; // trajTrackPairs_ now private
0457     const ConstTrajTrackPairCollection& trajtracks = eventInfo.trajTrackPairs();
0458 
0459     for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
0460          ++trajtrack) {
0461       m_counter_tracks++;
0462 
0463       const Trajectory* traj = (*trajtrack).first;
0464       const reco::Track* track = (*trajtrack).second;
0465 
0466       if (m_minTrackPt < track->pt() && track->pt() < m_maxTrackPt && m_minTrackP < track->p() &&
0467           track->p() < m_maxTrackP) {
0468         m_counter_trackmomentum++;
0469 
0470         if (fabs(track->dxy(eventInfo.beamSpot().position())) < m_maxDxy) {
0471           m_counter_trackdxy++;
0472           if (m_debug)
0473             std::cout << "JUST BEFORE muonResidualsFromTrack" << std::endl;
0474           MuonResidualsFromTrack muonResidualsFromTrack(builder,
0475                                                         magneticField,
0476                                                         globalGeometry,
0477                                                         muonDetIdAssociator,
0478                                                         prop,
0479                                                         traj,
0480                                                         track,
0481                                                         m_alignableNavigator,
0482                                                         1000.);
0483           if (m_debug)
0484             std::cout << "JUST AFTER muonResidualsFromTrack" << std::endl;
0485 
0486           if (m_debug)
0487             std::cout << "JUST BEFORE PROCESS" << std::endl;
0488           processMuonResidualsFromTrack(muonResidualsFromTrack);
0489           if (m_debug)
0490             std::cout << "JUST AFTER PROCESS" << std::endl;
0491         }
0492       }  // end if track p is within range
0493     }  // end if track pT is within range
0494     if (m_debug)
0495       std::cout << "JUST AFTER LOOP OVER trajTrackPairs" << std::endl;
0496 
0497   } else  // use muons
0498   {
0499     /*
0500            for (reco::MuonCollection::const_iterator muon = eventInfo.muonCollection_->begin();  muon != eventInfo.muonCollection_->end();  ++muon)
0501            {
0502            if ( !(muon->isTrackerMuon() && muon->innerTrack().isNonnull() ) ) continue;
0503 
0504            m_counter_tracks++;
0505 
0506            if (m_minTrackPt < muon->pt()  &&  muon->pt() < m_maxTrackPt && m_minTrackP < muon->p()  &&  muon->p() < m_maxTrackP)
0507            {
0508            m_counter_trackmomentum++;
0509 
0510            if (fabs(muon->innerTrack()->dxy(eventInfo.beamSpot_.position())) < m_maxDxy)
0511            {
0512            m_counter_trackdxy++;
0513 
0514         //std::cout<<"    *** will make MuonResidualsFromTrack ***"<<std::endl;
0515         MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), m_alignableNavigator, 100.);
0516         //std::cout<<"    *** have made MuonResidualsFromTrack ***"<<std::endl;
0517 
0518         //std::cout<<" trk eta="<<muon->eta()<<" ndof="<<muon->innerTrack()->ndof()<<" nchi2="<<muon->innerTrack()->normalizedChi2()
0519         //         <<" muresnchi2="<<muonResidualsFromTrack.normalizedChi2()<<" muresnhits="<<muonResidualsFromTrack.trackerNumHits()<<std::endl;
0520 
0521         processMuonResidualsFromTrack(muonResidualsFromTrack);
0522         } // end if track p is within range
0523         } // end if track pT is within range
0524         } // end loop over tracks
0525         */
0526   }
0527 }
0528 
0529 void MuonAlignmentFromReference::processMuonResidualsFromTrack(MuonResidualsFromTrack& mrft) {
0530   // std::cout << "minTrackerHits: " << mrft.trackerNumHits() << std::endl;
0531   if (mrft.trackerNumHits() >= m_minTrackerHits) {
0532     m_counter_trackerhits++;
0533     // std::cout << "mrft.normalizedChi2(): " << mrft.normalizedChi2() << std::endl;
0534 
0535     if (mrft.normalizedChi2() < m_maxTrackerRedChi2) {
0536       m_counter_trackerchi2++;
0537       if (m_allowTIDTEC || !mrft.contains_TIDTEC()) {
0538         m_counter_trackertidtec++;
0539 
0540         std::vector<DetId> chamberIds = mrft.chamberIds();
0541 
0542         if ((int)chamberIds.size() >= m_minNCrossedChambers) {
0543           m_counter_minchambers++;
0544 
0545           char charge = (mrft.getTrack()->charge() > 0 ? 1 : -1);
0546 
0547           for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end();
0548                ++chamberId) {
0549             if (chamberId->det() != DetId::Muon)
0550               continue;
0551             m_counter_totchambers++;
0552 
0553             // DT station 1,2,3
0554             if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT && DTChamberId(chamberId->rawId()).station() != 4) {
0555               MuonChamberResidual* dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
0556               MuonChamberResidual* dt2 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT2);
0557 
0558               m_counter_station123++;
0559               if (dt13 != nullptr && dt2 != nullptr) {
0560                 m_counter_station123valid++;
0561                 if (dt13->numHits() >= m_minDT13Hits) {
0562                   m_counter_station123dt13hits++;
0563                   if (dt2->numHits() >= m_minDT2Hits) {
0564                     m_counter_station123dt2hits++;
0565                     std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter =
0566                         m_fitters.find(dt13->chamberAlignable());
0567                     if (fitter != m_fitters.end()) {
0568                       m_counter_station123aligning++;
0569                       if (fabs(dt2->resslope()) < m_maxResSlopeY && (dt2->chi2() / double(dt2->ndof())) < 2.0) {
0570                         m_counter_resslopey++;
0571                         double* residdata = new double[MuonResiduals6DOFFitter::kNData];
0572                         residdata[MuonResiduals6DOFFitter::kResidX] = dt13->residual();
0573                         residdata[MuonResiduals6DOFFitter::kResidY] = dt2->residual();
0574                         residdata[MuonResiduals6DOFFitter::kResSlopeX] = dt13->resslope();
0575                         residdata[MuonResiduals6DOFFitter::kResSlopeY] = dt2->resslope();
0576                         residdata[MuonResiduals6DOFFitter::kPositionX] = dt13->trackx();
0577                         residdata[MuonResiduals6DOFFitter::kPositionY] = dt13->tracky();
0578                         residdata[MuonResiduals6DOFFitter::kAngleX] = dt13->trackdxdz();
0579                         residdata[MuonResiduals6DOFFitter::kAngleY] = dt13->trackdydz();
0580                         residdata[MuonResiduals6DOFFitter::kRedChi2] =
0581                             (dt13->chi2() + dt2->chi2()) / double(dt13->ndof() + dt2->ndof());
0582                         residdata[MuonResiduals6DOFFitter::kPz] = mrft.getTrack()->pz();
0583                         residdata[MuonResiduals6DOFFitter::kPt] = mrft.getTrack()->pt();
0584                         residdata[MuonResiduals6DOFFitter::kCharge] = mrft.getTrack()->charge();
0585                         residdata[MuonResiduals6DOFFitter::kStation] = DTChamberId(chamberId->rawId()).station();
0586                         residdata[MuonResiduals6DOFFitter::kWheel] = DTChamberId(chamberId->rawId()).wheel();
0587                         residdata[MuonResiduals6DOFFitter::kSector] = DTChamberId(chamberId->rawId()).sector();
0588                         residdata[MuonResiduals6DOFFitter::kChambW] = dt13->ChambW();
0589                         residdata[MuonResiduals6DOFFitter::kChambl] = dt13->Chambl();
0590 
0591                         if (m_debug) {
0592                           std::cout << "processMuonResidualsFromTrack 6DOF dt13->residual()  " << dt13->residual()
0593                                     << std::endl;
0594                           std::cout << "                                   dt2->residual()   " << dt2->residual()
0595                                     << std::endl;
0596                           std::cout << "                                   dt13->resslope()  " << dt13->resslope()
0597                                     << std::endl;
0598                           std::cout << "                                   dt2->resslope()   " << dt2->resslope()
0599                                     << std::endl;
0600                           std::cout << "                                   dt13->trackx()    " << dt13->trackx()
0601                                     << std::endl;
0602                           std::cout << "                                   dt13->tracky()    " << dt13->tracky()
0603                                     << std::endl;
0604                           std::cout << "                                   dt13->trackdxdz() " << dt13->trackdxdz()
0605                                     << std::endl;
0606                           std::cout << "                                   dt13->trackdydz() " << dt13->trackdydz()
0607                                     << std::endl;
0608                         }
0609 
0610                         fitter->second->fill(charge, residdata);
0611                         // the MuonResidualsFitter will delete the array when it is destroyed
0612                       }
0613                     }
0614                   }
0615                 }
0616               }
0617             }
0618 
0619             // DT 4th station
0620             else if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT &&
0621                      DTChamberId(chamberId->rawId()).station() == 4) {
0622               MuonChamberResidual* dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
0623 
0624               m_counter_station4++;
0625               if (dt13 != nullptr) {
0626                 m_counter_station4valid++;
0627                 if (dt13->numHits() >= m_minDT13Hits) {
0628                   m_counter_station4hits++;
0629 
0630                   std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter =
0631                       m_fitters.find(dt13->chamberAlignable());
0632                   if (fitter != m_fitters.end()) {
0633                     m_counter_station4aligning++;
0634 
0635                     double* residdata = new double[MuonResiduals5DOFFitter::kNData];
0636                     residdata[MuonResiduals5DOFFitter::kResid] = dt13->residual();
0637                     residdata[MuonResiduals5DOFFitter::kResSlope] = dt13->resslope();
0638                     residdata[MuonResiduals5DOFFitter::kPositionX] = dt13->trackx();
0639                     residdata[MuonResiduals5DOFFitter::kPositionY] = dt13->tracky();
0640                     residdata[MuonResiduals5DOFFitter::kAngleX] = dt13->trackdxdz();
0641                     residdata[MuonResiduals5DOFFitter::kAngleY] = dt13->trackdydz();
0642                     residdata[MuonResiduals5DOFFitter::kRedChi2] = dt13->chi2() / double(dt13->ndof());
0643                     residdata[MuonResiduals5DOFFitter::kPz] = mrft.getTrack()->pz();
0644                     residdata[MuonResiduals5DOFFitter::kPt] = mrft.getTrack()->pt();
0645                     residdata[MuonResiduals5DOFFitter::kCharge] = mrft.getTrack()->charge();
0646                     residdata[MuonResiduals5DOFFitter::kStation] = DTChamberId(chamberId->rawId()).station();
0647                     residdata[MuonResiduals5DOFFitter::kWheel] = DTChamberId(chamberId->rawId()).wheel();
0648                     residdata[MuonResiduals5DOFFitter::kSector] = DTChamberId(chamberId->rawId()).sector();
0649                     residdata[MuonResiduals5DOFFitter::kChambW] = dt13->ChambW();
0650                     residdata[MuonResiduals5DOFFitter::kChambl] = dt13->Chambl();
0651 
0652                     if (m_debug) {
0653                       std::cout << "processMuonResidualsFromTrack 5DOF dt13->residual()  " << dt13->residual()
0654                                 << std::endl;
0655                       std::cout << "                                   dt13->resslope()  " << dt13->resslope()
0656                                 << std::endl;
0657                       std::cout << "                                   dt13->trackx()    " << dt13->trackx()
0658                                 << std::endl;
0659                       std::cout << "                                   dt13->tracky()    " << dt13->tracky()
0660                                 << std::endl;
0661                       std::cout << "                                   dt13->trackdxdz() " << dt13->trackdxdz()
0662                                 << std::endl;
0663                       std::cout << "                                   dt13->trackdydz() " << dt13->trackdydz()
0664                                 << std::endl;
0665                     }
0666 
0667                     fitter->second->fill(charge, residdata);
0668                     // the MuonResidualsFitter will delete the array when it is destroyed
0669                   }
0670                 }
0671               }
0672             }  // end DT 4th station
0673 
0674             // CSC
0675             else if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC) {
0676               MuonChamberResidual* csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
0677               m_counter_csc++;
0678               if (csc != nullptr) {
0679                 m_counter_cscvalid++;
0680                 if (csc->numHits() >= m_minCSCHits) {
0681                   m_counter_cschits++;
0682                   Alignable* ali = csc->chamberAlignable();
0683 
0684                   CSCDetId id(ali->geomDetId().rawId());
0685                   if (m_combineME11 && id.station() == 1 && id.ring() == 4)
0686                     ali = m_me11map[ali];
0687 
0688                   std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(ali);
0689                   if (fitter != m_fitters.end()) {
0690                     m_counter_cscaligning++;
0691                     double* residdata = new double[MuonResiduals6DOFrphiFitter::kNData];
0692                     residdata[MuonResiduals6DOFrphiFitter::kResid] = csc->residual();
0693                     residdata[MuonResiduals6DOFrphiFitter::kResSlope] = csc->resslope();
0694                     residdata[MuonResiduals6DOFrphiFitter::kPositionX] = csc->trackx();
0695                     residdata[MuonResiduals6DOFrphiFitter::kPositionY] = csc->tracky();
0696                     residdata[MuonResiduals6DOFrphiFitter::kAngleX] = csc->trackdxdz();
0697                     residdata[MuonResiduals6DOFrphiFitter::kAngleY] = csc->trackdydz();
0698                     residdata[MuonResiduals6DOFrphiFitter::kRedChi2] = csc->chi2() / double(csc->ndof());
0699                     residdata[MuonResiduals6DOFrphiFitter::kPz] = mrft.getTrack()->pz();
0700                     residdata[MuonResiduals6DOFrphiFitter::kPt] = mrft.getTrack()->pt();
0701                     residdata[MuonResiduals6DOFrphiFitter::kCharge] = mrft.getTrack()->charge();
0702 
0703                     if (m_debug) {
0704                       std::cout << "processMuonResidualsFromTrack 6DOFrphi csc->residual()  " << csc->residual()
0705                                 << std::endl;
0706                       std::cout << "                                       csc->resslope()  " << csc->resslope()
0707                                 << std::endl;
0708                       std::cout << "                                       csc->trackx()    " << csc->trackx()
0709                                 << std::endl;
0710                       std::cout << "                                       csc->tracky()    " << csc->tracky()
0711                                 << std::endl;
0712                       std::cout << "                                       csc->trackdxdz() " << csc->trackdxdz()
0713                                 << std::endl;
0714                       std::cout << "                                       csc->trackdydz() " << csc->trackdydz()
0715                                 << std::endl;
0716                     }
0717 
0718                     fitter->second->fill(charge, residdata);
0719                     // the MuonResidualsFitter will delete the array when it is destroyed
0720                   }
0721                 }
0722               }
0723             }  // end CSC
0724 
0725             else if (m_doDT && m_doCSC)
0726               assert(false);
0727 
0728           }  // end loop over chamberIds
0729         }  // # crossed muon chambers ok
0730       }  // endcap tracker ok
0731     }  // chi2 ok
0732   }  // trackerNumHits ok
0733 }
0734 
0735 void MuonAlignmentFromReference::terminate(const edm::EventSetup& iSetup) {
0736   bool m_debug = false;
0737 
0738   // one-time print-out
0739   std::cout << "Counters:" << std::endl
0740             << "COUNT{ events: " << m_counter_events << " }" << std::endl
0741             << "COUNT{  tracks: " << m_counter_tracks << " }" << std::endl
0742             << "COUNT{   trackppt: " << m_counter_trackmomentum << " }" << std::endl
0743             << "COUNT{    trackdxy: " << m_counter_trackdxy << " }" << std::endl
0744             << "COUNT{     trackerhits: " << m_counter_trackerhits << " }" << std::endl
0745             << "COUNT{      trackerchi2: " << m_counter_trackerchi2 << " }" << std::endl
0746             << "COUNT{       trackertidtec: " << m_counter_trackertidtec << " }" << std::endl
0747             << "COUNT{        minnchambers: " << m_counter_minchambers << " }" << std::endl
0748             << "COUNT{          totchambers: " << m_counter_totchambers << " }" << std::endl
0749             << "COUNT{            station123:             " << m_counter_station123 << " }" << std::endl
0750             << "COUNT{             station123valid:       " << m_counter_station123valid << " }" << std::endl
0751             << "COUNT{              station123dt13hits:   " << m_counter_station123dt13hits << " }" << std::endl
0752             << "COUNT{               station123dt2hits:   " << m_counter_station123dt2hits << " }" << std::endl
0753             << "COUNT{                station123aligning: " << m_counter_station123aligning << " }" << std::endl
0754             << "COUNT{                 resslopey: " << m_counter_resslopey << " }" << std::endl
0755             << "COUNT{            station4:            " << m_counter_station4 << " }" << std::endl
0756             << "COUNT{             station4valid:      " << m_counter_station4valid << " }" << std::endl
0757             << "COUNT{              station4hits:      " << m_counter_station4hits << " }" << std::endl
0758             << "COUNT{               station4aligning: " << m_counter_station4aligning << " }" << std::endl
0759             << "COUNT{            csc:            " << m_counter_csc << " }" << std::endl
0760             << "COUNT{             cscvalid:      " << m_counter_cscvalid << " }" << std::endl
0761             << "COUNT{              cschits:      " << m_counter_cschits << " }" << std::endl
0762             << "COUNT{               cscaligning: " << m_counter_cscaligning << " }" << std::endl
0763             << "That's all!" << std::endl;
0764 
0765   TStopwatch stop_watch;
0766 
0767   // collect temporary files
0768   if (!m_readTemporaryFiles.empty()) {
0769     stop_watch.Start();
0770     readTmpFiles();
0771     if (m_debug)
0772       std::cout << "readTmpFiles took " << stop_watch.CpuTime() << " sec" << std::endl;
0773     stop_watch.Stop();
0774   }
0775 
0776   // select residuals peaks and discard tails if peakNSigma>0 (only while doing alignment)
0777   if (m_peakNSigma > 0. && m_doAlignment) {
0778     stop_watch.Start();
0779     selectResidualsPeaks();
0780     if (m_debug)
0781       std::cout << "selectResidualsPeaks took " << stop_watch.CpuTime() << " sec" << std::endl;
0782     stop_watch.Stop();
0783   }
0784 
0785   if (m_BFieldCorrection > 0 && m_doAlignment) {
0786     stop_watch.Start();
0787     correctBField();
0788     if (m_debug)
0789       std::cout << "correctBField took " << stop_watch.CpuTime() << " sec" << std::endl;
0790     stop_watch.Stop();
0791   }
0792 
0793   if (m_doAlignment && !m_doCSC)  // for now apply fiducial cuts to DT only
0794   {
0795     stop_watch.Start();
0796     fiducialCuts();
0797     if (m_debug)
0798       std::cout << "fiducialCuts took " << stop_watch.CpuTime() << " sec" << std::endl;
0799     stop_watch.Stop();
0800   }
0801 
0802   // optionally, create an nutuple for easy debugging
0803   if (m_createNtuple) {
0804     stop_watch.Start();
0805     fillNtuple();
0806     if (m_debug)
0807       std::cout << "fillNtuple took " << stop_watch.CpuTime() << " sec" << std::endl;
0808     stop_watch.Stop();
0809   }
0810 
0811   if (m_doAlignment) {
0812     stop_watch.Start();
0813     eraseNotSelectedResiduals();
0814     if (m_debug)
0815       std::cout << "eraseNotSelectedResiduals took " << stop_watch.CpuTime() << " sec" << std::endl;
0816     stop_watch.Stop();
0817   }
0818 
0819   // fit and align (time-consuming, so the user can turn it off if in a residuals-gathering job)
0820   if (m_doAlignment) {
0821     stop_watch.Start();
0822     fitAndAlign();
0823     if (m_debug)
0824       std::cout << "fitAndAlign took " << stop_watch.CpuTime() << " sec" << std::endl;
0825     stop_watch.Stop();
0826   }
0827 
0828   // write out the pseudontuples for a later job to collect
0829   if (m_writeTemporaryFile != std::string(""))
0830     writeTmpFiles();
0831   if (m_debug)
0832     std::cout << "end: MuonAlignmentFromReference::terminate()" << std::endl;
0833 }
0834 
0835 void MuonAlignmentFromReference::fitAndAlign() {
0836   bool m_debug = false;
0837 
0838   edm::Service<TFileService> tfileService;
0839   TFileDirectory rootDirectory(tfileService->mkdir("MuonAlignmentFromReference"));
0840 
0841   std::ofstream report;
0842   bool writeReport = (m_reportFileName != std::string(""));
0843   if (writeReport) {
0844     report.open(m_reportFileName.c_str());
0845     report << "nan = None;  NAN = None" << std::endl;
0846     report << "nan = 0" << std::endl;
0847     report << "reports = []" << std::endl;
0848     report << "class ValErr:" << std::endl
0849            << "    def __init__(self, value, error, antisym):" << std::endl
0850            << "        self.value, self.error, self.antisym = value, error, antisym" << std::endl
0851            << "" << std::endl
0852            << "    def __repr__(self):" << std::endl
0853            << "        if self.antisym == 0.:" << std::endl
0854            << "            return \"%g +- %g\" % (self.value, self.error)" << std::endl
0855            << "        else:" << std::endl
0856            << "            return \"%g +- %g ~ %g\" % (self.value, self.error, self.antisym)" << std::endl
0857            << "" << std::endl
0858            << "class Report:" << std::endl
0859            << "    def __init__(self, chamberId, postal_address, name):" << std::endl
0860            << "        self.chamberId, self.postal_address, self.name = chamberId, postal_address, name" << std::endl
0861            << "        self.status = \"NOFIT\"" << std::endl
0862            << "        self.fittype = None" << std::endl
0863            << "" << std::endl
0864            << "    def add_parameters(self, deltax, deltay, deltaz, deltaphix, deltaphiy, deltaphiz, loglikelihood, "
0865               "numsegments, sumofweights, redchi2):"
0866            << std::endl
0867            << "        self.status = \"PASS\"" << std::endl
0868            << "        self.deltax, self.deltay, self.deltaz, self.deltaphix, self.deltaphiy, self.deltaphiz = deltax, "
0869               "deltay, deltaz, deltaphix, deltaphiy, deltaphiz"
0870            << std::endl
0871            << "        self.loglikelihood, self.numsegments, self.sumofweights, self.redchi2 = loglikelihood, "
0872               "numsegments, sumofweights, redchi2"
0873            << std::endl
0874            << "" << std::endl
0875            << "    def add_stats(self, median_x, median_y, median_dxdz, median_dydz, mean30_x, mean30_y, mean20_dxdz, "
0876               "mean50_dydz, mean15_x, mean15_y, mean10_dxdz, mean25_dydz, wmean30_x, wmean30_y, wmean20_dxdz, "
0877               "wmean50_dydz, wmean15_x, wmean15_y, wmean10_dxdz, wmean25_dydz, stdev30_x, stdev30_y, stdev20_dxdz, "
0878               "stdev50_dydz, stdev15_x, stdev15_y, stdev10_dxdz, stdev25_dydz):"
0879            << std::endl
0880            << "        self.median_x, self.median_y, self.median_dxdz, self.median_dydz, self.mean30_x, self.mean30_y, "
0881               "self.mean20_dxdz, self.mean50_dydz, self.mean15_x, self.mean15_y, self.mean10_dxdz, self.mean25_dydz, "
0882               "self.wmean30_x, self.wmean30_y, self.wmean20_dxdz, self.wmean50_dydz, self.wmean15_x, self.wmean15_y, "
0883               "self.wmean10_dxdz, self.wmean25_dydz, self.stdev30_x, self.stdev30_y, self.stdev20_dxdz, "
0884               "self.stdev50_dydz, self.stdev15_x, self.stdev15_y, self.stdev10_dxdz, self.stdev25_dydz = median_x, "
0885               "median_y, median_dxdz, median_dydz, mean30_x, mean30_y, mean20_dxdz, mean50_dydz, mean15_x, mean15_y, "
0886               "mean10_dxdz, mean25_dydz, wmean30_x, wmean30_y, wmean20_dxdz, wmean50_dydz, wmean15_x, wmean15_y, "
0887               "wmean10_dxdz, wmean25_dydz, stdev30_x, stdev30_y, stdev20_dxdz, stdev50_dydz, stdev15_x, stdev15_y, "
0888               "stdev10_dxdz, stdev25_dydz"
0889            << std::endl
0890            << "" << std::endl
0891            << "    def __repr__(self):" << std::endl
0892            << "        return \"<Report %s %s %s>\" % (self.postal_address[0], \" \".join(map(str, "
0893               "self.postal_address[1:])), self.status)"
0894            << std::endl
0895            << std::endl;
0896   }
0897 
0898   if (m_debug)
0899     std::cout << "***** just after report.open" << std::endl;
0900 
0901   for (const auto& ali : m_alignables) {
0902     if (m_debug)
0903       std::cout << "***** Start loop over alignables" << std::endl;
0904 
0905     std::vector<bool> selector = ali->alignmentParameters()->selector();
0906     bool align_x = selector[0];
0907     bool align_y = selector[1];
0908     bool align_z = selector[2];
0909     bool align_phix = selector[3];
0910     bool align_phiy = selector[4];
0911     bool align_phiz = selector[5];
0912     int numParams = ((align_x ? 1 : 0) + (align_y ? 1 : 0) + (align_z ? 1 : 0) + (align_phix ? 1 : 0) +
0913                      (align_phiy ? 1 : 0) + (align_phiz ? 1 : 0));
0914 
0915     // map from 0-5 to the index of params, above
0916     std::vector<int> paramIndex;
0917     int paramIndex_counter = -1;
0918     if (align_x)
0919       paramIndex_counter++;
0920     paramIndex.push_back(paramIndex_counter);
0921     if (align_y)
0922       paramIndex_counter++;
0923     paramIndex.push_back(paramIndex_counter);
0924     if (align_z)
0925       paramIndex_counter++;
0926     paramIndex.push_back(paramIndex_counter);
0927     if (align_phix)
0928       paramIndex_counter++;
0929     paramIndex.push_back(paramIndex_counter);
0930     if (align_phiy)
0931       paramIndex_counter++;
0932     paramIndex.push_back(paramIndex_counter);
0933     if (align_phiz)
0934       paramIndex_counter++;
0935     paramIndex.push_back(paramIndex_counter);
0936 
0937     DetId id = ali->geomDetId();
0938 
0939     auto thisali = ali;
0940     if (m_combineME11 && id.subdetId() == MuonSubdetId::CSC) {
0941       CSCDetId cscid(id.rawId());
0942       if (cscid.station() == 1 && cscid.ring() == 4)
0943         thisali = m_me11map[ali];
0944     }
0945 
0946     if (m_debug)
0947       std::cout << "***** loop over alignables 1" << std::endl;
0948 
0949     char cname[40];
0950     char wheel_label[][2] = {"A", "B", "C", "D", "E"};
0951 
0952     if (id.subdetId() == MuonSubdetId::DT) {
0953       DTChamberId chamberId(id.rawId());
0954 
0955       //if ( ! ( (chamberId.station()==1&&chamberId.wheel()==0) || (chamberId.station()==4&&chamberId.wheel()==2) ) ) continue;
0956 
0957       sprintf(cname, "MBwh%sst%dsec%02d", wheel_label[chamberId.wheel() + 2], chamberId.station(), chamberId.sector());
0958       if (writeReport) {
0959         report << "reports.append(Report(" << id.rawId() << ", (\"DT\", " << chamberId.wheel() << ", "
0960                << chamberId.station() << ", " << chamberId.sector() << "), \"" << cname << "\"))" << std::endl;
0961       }
0962     } else if (id.subdetId() == MuonSubdetId::CSC) {
0963       CSCDetId chamberId(id.rawId());
0964       sprintf(cname,
0965               "ME%s%d%d_%02d",
0966               (chamberId.endcap() == 1 ? "p" : "m"),
0967               chamberId.station(),
0968               chamberId.ring(),
0969               chamberId.chamber());
0970 
0971       //if ( chamberId.chamber()>6 || chamberId.endcap()==2 || ! ( (chamberId.station()==2&&chamberId.ring()==1) || (chamberId.station()==3&&chamberId.ring()==2) ) ) continue;
0972 
0973       if (writeReport) {
0974         report << "reports.append(Report(" << id.rawId() << ", (\"CSC\", " << chamberId.endcap() << ", "
0975                << chamberId.station() << ", " << chamberId.ring() << ", " << chamberId.chamber() << "), \"" << cname
0976                << "\"))" << std::endl;
0977       }
0978     }
0979 
0980     if (m_debug)
0981       std::cout << "***** loop over alignables 2" << std::endl;
0982 
0983     //if(! ( strcmp(cname,"MBwhCst3sec12")==0 || strcmp(cname,"MBwhCst3sec06")==0)) continue;
0984 
0985     std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(thisali);
0986 
0987     if (m_debug)
0988       std::cout << "***** loop over alignables 3" << std::endl;
0989 
0990     if (fitter != m_fitters.end()) {
0991       //if (fitter->second->type() != MuonResidualsFitter::k6DOFrphi) continue;
0992 
0993       TStopwatch stop_watch;
0994       stop_watch.Start();
0995 
0996       // MINUIT is verbose in std::cout anyway
0997       if (m_debug)
0998         std::cout << "============================================================================================="
0999                   << std::endl;
1000       if (m_debug)
1001         std::cout << "Fitting " << cname << std::endl;
1002 
1003       if (writeReport) {
1004         report << "reports[-1].posNum = " << fitter->second->numResidualsPos() << std::endl;
1005         report << "reports[-1].negNum = " << fitter->second->numResidualsNeg() << std::endl;
1006       }
1007 
1008       if (fitter->second->type() == MuonResidualsFitter::k5DOF) {
1009         if (!align_x)
1010           fitter->second->fix(MuonResiduals5DOFFitter::kAlignX);
1011         if (!align_z)
1012           fitter->second->fix(MuonResiduals5DOFFitter::kAlignZ);
1013         if (!align_phix)
1014           fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiX);
1015         if (!align_phiy)
1016           fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiY);
1017         if (!align_phiz)
1018           fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiZ);
1019       } else if (fitter->second->type() == MuonResidualsFitter::k6DOF) {
1020         if (!align_x)
1021           fitter->second->fix(MuonResiduals6DOFFitter::kAlignX);
1022         if (!align_y)
1023           fitter->second->fix(MuonResiduals6DOFFitter::kAlignY);
1024         if (!align_z)
1025           fitter->second->fix(MuonResiduals6DOFFitter::kAlignZ);
1026         if (!align_phix)
1027           fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiX);
1028         if (!align_phiy)
1029           fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiY);
1030         if (!align_phiz)
1031           fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiZ);
1032       } else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi) {
1033         if (!align_x)
1034           fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignX);
1035         if (!align_y)
1036           fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignY);
1037         if (!align_z)
1038           fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignZ);
1039         if (!align_phix)
1040           fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1041         if (!align_phiy)
1042           fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1043         if (!align_phiz)
1044           fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1045       } else
1046         assert(false);
1047 
1048       if (m_debug)
1049         std::cout << "***** loop over alignables 4" << std::endl;
1050 
1051       AlgebraicVector params(numParams);
1052       AlgebraicSymMatrix cov(numParams);
1053 
1054       if (fitter->second->numsegments() >= m_minAlignmentHits) {
1055         if (m_debug)
1056           std::cout << "***** loop over alignables 5" << std::endl;
1057 
1058         bool successful_fit = fitter->second->fit(thisali);
1059 
1060         if (m_debug)
1061           std::cout << "***** loop over alignables 6 " << fitter->second->type() << std::endl;
1062 
1063         double loglikelihood = fitter->second->loglikelihood();
1064         double numsegments = fitter->second->numsegments();
1065         double sumofweights = fitter->second->sumofweights();
1066         double redchi2 = fitter->second->plot(cname, &rootDirectory, thisali);
1067 
1068         if (fitter->second->type() == MuonResidualsFitter::k5DOF) {
1069           if (m_debug)
1070             std::cout << "***** loop over alignables k5DOF" << std::endl;
1071 
1072           double deltax_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignX);
1073           double deltax_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignX);
1074           double deltax_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignX);
1075 
1076           double deltaz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignZ);
1077           double deltaz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignZ);
1078           double deltaz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignZ);
1079 
1080           double deltaphix_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiX);
1081           double deltaphix_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiX);
1082           double deltaphix_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiX);
1083 
1084           double deltaphiy_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiY);
1085           double deltaphiy_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiY);
1086           double deltaphiy_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiY);
1087 
1088           double deltaphiz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiZ);
1089           double deltaphiz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiZ);
1090           double deltaphiz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiZ);
1091 
1092           double sigmaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidSigma);
1093           double sigmaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidSigma);
1094           double sigmaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidSigma);
1095 
1096           double sigmaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeSigma);
1097           double sigmaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeSigma);
1098           double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeSigma);
1099 
1100           double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error,
1101               gammaresslope_antisym;
1102           gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error =
1103               gammaresslope_antisym = 0.;
1104 
1105           if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1106               fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1107               fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1108             gammaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidGamma);
1109             gammaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidGamma);
1110             gammaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidGamma);
1111 
1112             gammaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeGamma);
1113             gammaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeGamma);
1114             gammaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeGamma);
1115           }
1116 
1117           if (writeReport) {
1118             report << "reports[-1].fittype = \"5DOF\"" << std::endl;
1119             report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", "
1120                    << deltax_antisym << "), \\" << std::endl
1121                    << "                           None, \\" << std::endl
1122                    << "                           ValErr(" << deltaz_value << ", " << deltaz_error << ", "
1123                    << deltaz_antisym << "), \\" << std::endl
1124                    << "                           ValErr(" << deltaphix_value << ", " << deltaphix_error << ", "
1125                    << deltaphix_antisym << "), \\" << std::endl
1126                    << "                           ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", "
1127                    << deltaphiy_antisym << "), \\" << std::endl
1128                    << "                           ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", "
1129                    << deltaphiz_antisym << "), \\" << std::endl
1130                    << "                           " << loglikelihood << ", " << numsegments << ", " << sumofweights
1131                    << ", " << redchi2 << ")" << std::endl;
1132             report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", "
1133                    << sigmaresid_antisym << ")" << std::endl;
1134             report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error
1135                    << ", " << sigmaresslope_antisym << ")" << std::endl;
1136             if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1137                 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1138                 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1139               report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", "
1140                      << gammaresid_antisym << ")" << std::endl;
1141               report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error
1142                      << ", " << gammaresslope_antisym << ")" << std::endl;
1143             }
1144 
1145             report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals5DOFFitter::kResid) << ", "
1146                    << "None, " << fitter->second->median(MuonResiduals5DOFFitter::kResSlope) << ", "
1147                    << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 30.) << ", "
1148                    << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 20.) << ", "
1149                    << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 15.) << ", "
1150                    << "None, " << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 10.) << ", "
1151                    << "None, "
1152                    << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 30.)
1153                    << ", "
1154                    << "None, "
1155                    << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 20.)
1156                    << ", "
1157                    << "None, "
1158                    << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 15.)
1159                    << ", "
1160                    << "None, "
1161                    << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 10.)
1162                    << ", "
1163                    << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 30.) << ", "
1164                    << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 20.) << ", "
1165                    << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 15.) << ", "
1166                    << "None, " << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 10.) << ", "
1167                    << "None)" << std::endl;
1168 
1169             std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
1170             namesimple_x << cname << "_simple_x";
1171             namesimple_dxdz << cname << "_simple_dxdz";
1172             nameweighted_x << cname << "_weighted_x";
1173             nameweighted_dxdz << cname << "_weighted_dxdz";
1174 
1175             fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals5DOFFitter::kResid, 10.);
1176             fitter->second->plotsimple(
1177                 namesimple_dxdz.str(), &rootDirectory, MuonResiduals5DOFFitter::kResSlope, 1000.);
1178 
1179             fitter->second->plotweighted(nameweighted_x.str(),
1180                                          &rootDirectory,
1181                                          MuonResiduals5DOFFitter::kResid,
1182                                          MuonResiduals5DOFFitter::kRedChi2,
1183                                          10.);
1184             fitter->second->plotweighted(nameweighted_dxdz.str(),
1185                                          &rootDirectory,
1186                                          MuonResiduals5DOFFitter::kResSlope,
1187                                          MuonResiduals5DOFFitter::kRedChi2,
1188                                          1000.);
1189           }
1190 
1191           if (successful_fit) {
1192             if (align_x)
1193               params[paramIndex[0]] = deltax_value;
1194             if (align_z)
1195               params[paramIndex[2]] = deltaz_value;
1196             if (align_phix)
1197               params[paramIndex[3]] = deltaphix_value;
1198             if (align_phiy)
1199               params[paramIndex[4]] = deltaphiy_value;
1200             if (align_phiz)
1201               params[paramIndex[5]] = deltaphiz_value;
1202           }
1203         }  // end if 5DOF
1204 
1205         else if (fitter->second->type() == MuonResidualsFitter::k6DOF) {
1206           if (m_debug)
1207             std::cout << "***** loop over alignables k6DOF" << std::endl;
1208 
1209           double deltax_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignX);
1210           double deltax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignX);
1211           double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignX);
1212 
1213           double deltay_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignY);
1214           double deltay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignY);
1215           double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignY);
1216 
1217           double deltaz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignZ);
1218           double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignZ);
1219           double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignZ);
1220 
1221           double deltaphix_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiX);
1222           double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiX);
1223           double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiX);
1224 
1225           double deltaphiy_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiY);
1226           double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiY);
1227           double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiY);
1228 
1229           double deltaphiz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiZ);
1230           double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiZ);
1231           double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiZ);
1232 
1233           double sigmax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXSigma);
1234           double sigmax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXSigma);
1235           double sigmax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXSigma);
1236 
1237           double sigmay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYSigma);
1238           double sigmay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYSigma);
1239           double sigmay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYSigma);
1240 
1241           double sigmadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXSigma);
1242           double sigmadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXSigma);
1243           double sigmadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXSigma);
1244 
1245           double sigmadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYSigma);
1246           double sigmadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYSigma);
1247           double sigmadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYSigma);
1248 
1249           double gammax_value, gammax_error, gammax_antisym, gammay_value, gammay_error, gammay_antisym,
1250               gammadxdz_value, gammadxdz_error, gammadxdz_antisym, gammadydz_value, gammadydz_error, gammadydz_antisym;
1251           gammax_value = gammax_error = gammax_antisym = gammay_value = gammay_error = gammay_antisym =
1252               gammadxdz_value = gammadxdz_error = gammadxdz_antisym = gammadydz_value = gammadydz_error =
1253                   gammadydz_antisym = 0.;
1254           if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1255               fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1256               fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1257             gammax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXGamma);
1258             gammax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXGamma);
1259             gammax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXGamma);
1260 
1261             gammay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYGamma);
1262             gammay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYGamma);
1263             gammay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYGamma);
1264 
1265             gammadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXGamma);
1266             gammadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXGamma);
1267             gammadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXGamma);
1268 
1269             gammadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYGamma);
1270             gammadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYGamma);
1271             gammadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYGamma);
1272           }
1273 
1274           if (writeReport) {
1275             report << "reports[-1].fittype = \"6DOF\"" << std::endl;
1276             report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", "
1277                    << deltax_antisym << "), \\" << std::endl
1278                    << "                           ValErr(" << deltay_value << ", " << deltay_error << ", "
1279                    << deltay_antisym << "), \\" << std::endl
1280                    << "                           ValErr(" << deltaz_value << ", " << deltaz_error << ", "
1281                    << deltaz_antisym << "), \\" << std::endl
1282                    << "                           ValErr(" << deltaphix_value << ", " << deltaphix_error << ", "
1283                    << deltaphix_antisym << "), \\" << std::endl
1284                    << "                           ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", "
1285                    << deltaphiy_antisym << "), \\" << std::endl
1286                    << "                           ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", "
1287                    << deltaphiz_antisym << "), \\" << std::endl
1288                    << "                           " << loglikelihood << ", " << numsegments << ", " << sumofweights
1289                    << ", " << redchi2 << ")" << std::endl;
1290             report << "reports[-1].sigmax = ValErr(" << sigmax_value << ", " << sigmax_error << ", " << sigmax_antisym
1291                    << ")" << std::endl;
1292             report << "reports[-1].sigmay = ValErr(" << sigmay_value << ", " << sigmay_error << ", " << sigmay_antisym
1293                    << ")" << std::endl;
1294             report << "reports[-1].sigmadxdz = ValErr(" << sigmadxdz_value << ", " << sigmadxdz_error << ", "
1295                    << sigmadxdz_antisym << ")" << std::endl;
1296             report << "reports[-1].sigmadydz = ValErr(" << sigmadydz_value << ", " << sigmadydz_error << ", "
1297                    << sigmadydz_antisym << ")" << std::endl;
1298             if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1299                 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1300                 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1301               report << "reports[-1].gammax = ValErr(" << gammax_value << ", " << gammax_error << ", " << gammax_antisym
1302                      << ")" << std::endl;
1303               report << "reports[-1].gammay = ValErr(" << gammay_value << ", " << gammay_error << ", " << gammay_antisym
1304                      << ")" << std::endl;
1305               report << "reports[-1].gammadxdz = ValErr(" << gammadxdz_value << ", " << gammadxdz_error << ", "
1306                      << gammadxdz_antisym << ")" << std::endl;
1307               report << "reports[-1].gammadydz = ValErr(" << gammadydz_value << ", " << gammadydz_error << ", "
1308                      << gammadydz_antisym << ")" << std::endl;
1309             }
1310 
1311             report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals6DOFFitter::kResidX) << ", "
1312                    << fitter->second->median(MuonResiduals6DOFFitter::kResidY) << ", "
1313                    << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeX) << ", "
1314                    << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeY) << ", "
1315                    << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
1316                    << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
1317                    << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
1318                    << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
1319                    << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
1320                    << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
1321                    << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
1322                    << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ", "
1323                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 30.)
1324                    << ", "
1325                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 30.)
1326                    << ", "
1327                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 20.)
1328                    << ", "
1329                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 50.)
1330                    << ", "
1331                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 15.)
1332                    << ", "
1333                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 15.)
1334                    << ", "
1335                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 10.)
1336                    << ", "
1337                    << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 25.)
1338                    << ", " << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
1339                    << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
1340                    << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
1341                    << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
1342                    << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
1343                    << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
1344                    << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
1345                    << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ")" << std::endl;
1346 
1347             std::stringstream namesimple_x, namesimple_y, namesimple_dxdz, namesimple_dydz, nameweighted_x,
1348                 nameweighted_y, nameweighted_dxdz, nameweighted_dydz;
1349             namesimple_x << cname << "_simple_x";
1350             namesimple_y << cname << "_simple_y";
1351             namesimple_dxdz << cname << "_simple_dxdz";
1352             namesimple_dydz << cname << "_simple_dydz";
1353             nameweighted_x << cname << "_weighted_x";
1354             nameweighted_y << cname << "_weighted_y";
1355             nameweighted_dxdz << cname << "_weighted_dxdz";
1356             nameweighted_dydz << cname << "_weighted_dydz";
1357 
1358             fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidX, 10.);
1359             fitter->second->plotsimple(namesimple_y.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidY, 10.);
1360             fitter->second->plotsimple(
1361                 namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeX, 1000.);
1362             fitter->second->plotsimple(
1363                 namesimple_dydz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeY, 1000.);
1364 
1365             fitter->second->plotweighted(nameweighted_x.str(),
1366                                          &rootDirectory,
1367                                          MuonResiduals6DOFFitter::kResidX,
1368                                          MuonResiduals6DOFFitter::kRedChi2,
1369                                          10.);
1370             fitter->second->plotweighted(nameweighted_y.str(),
1371                                          &rootDirectory,
1372                                          MuonResiduals6DOFFitter::kResidY,
1373                                          MuonResiduals6DOFFitter::kRedChi2,
1374                                          10.);
1375             fitter->second->plotweighted(nameweighted_dxdz.str(),
1376                                          &rootDirectory,
1377                                          MuonResiduals6DOFFitter::kResSlopeX,
1378                                          MuonResiduals6DOFFitter::kRedChi2,
1379                                          1000.);
1380             fitter->second->plotweighted(nameweighted_dydz.str(),
1381                                          &rootDirectory,
1382                                          MuonResiduals6DOFFitter::kResSlopeY,
1383                                          MuonResiduals6DOFFitter::kRedChi2,
1384                                          1000.);
1385           }
1386 
1387           if (successful_fit) {
1388             if (align_x)
1389               params[paramIndex[0]] = deltax_value;
1390             if (align_y)
1391               params[paramIndex[1]] = deltay_value;
1392             if (align_z)
1393               params[paramIndex[2]] = deltaz_value;
1394             if (align_phix)
1395               params[paramIndex[3]] = deltaphix_value;
1396             if (align_phiy)
1397               params[paramIndex[4]] = deltaphiy_value;
1398             if (align_phiz)
1399               params[paramIndex[5]] = deltaphiz_value;
1400           }
1401         }  // end if 6DOF
1402 
1403         else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi) {
1404           if (m_debug)
1405             std::cout << "***** loop over alignables k6DOFrphi" << std::endl;
1406 
1407           double deltax_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignX);
1408           double deltax_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignX);
1409           double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignX);
1410 
1411           double deltay_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignY);
1412           double deltay_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignY);
1413           double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignY);
1414 
1415           double deltaz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignZ);
1416           double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignZ);
1417           double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignZ);
1418 
1419           double deltaphix_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1420           double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1421           double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiX);
1422 
1423           double deltaphiy_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1424           double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1425           double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiY);
1426 
1427           double deltaphiz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1428           double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1429           double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
1430 
1431           double sigmaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidSigma);
1432           double sigmaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidSigma);
1433           double sigmaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidSigma);
1434 
1435           double sigmaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
1436           double sigmaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
1437           double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
1438 
1439           double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error,
1440               gammaresslope_antisym;
1441           gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error =
1442               gammaresslope_antisym = 0.;
1443           if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1444               fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1445               fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1446             gammaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidGamma);
1447             gammaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidGamma);
1448             gammaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidGamma);
1449 
1450             gammaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
1451             gammaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
1452             gammaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
1453           }
1454 
1455           if (writeReport) {
1456             report << "reports[-1].fittype = \"6DOFrphi\"" << std::endl;
1457             report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", "
1458                    << deltax_antisym << "), \\" << std::endl
1459                    << "                           ValErr(" << deltay_value << ", " << deltay_error << ", "
1460                    << deltay_antisym << "), \\" << std::endl
1461                    << "                           ValErr(" << deltaz_value << ", " << deltaz_error << ", "
1462                    << deltaz_antisym << "), \\" << std::endl
1463                    << "                           ValErr(" << deltaphix_value << ", " << deltaphix_error << ", "
1464                    << deltaphix_antisym << "), \\" << std::endl
1465                    << "                           ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", "
1466                    << deltaphiy_antisym << "), \\" << std::endl
1467                    << "                           ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", "
1468                    << deltaphiz_antisym << "), \\" << std::endl
1469                    << "                           " << loglikelihood << ", " << numsegments << ", " << sumofweights
1470                    << ", " << redchi2 << ")" << std::endl;
1471             report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", "
1472                    << sigmaresid_antisym << ")" << std::endl;
1473             report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error
1474                    << ", " << sigmaresslope_antisym << ")" << std::endl;
1475             if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
1476                 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
1477                 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
1478               report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", "
1479                      << gammaresid_antisym << ")" << std::endl;
1480               report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error
1481                      << ", " << gammaresslope_antisym << ")" << std::endl;
1482             }
1483 
1484             report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals6DOFrphiFitter::kResid) << ", "
1485                    << "None, " << fitter->second->median(MuonResiduals6DOFrphiFitter::kResSlope) << ", "
1486                    << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", "
1487                    << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", "
1488                    << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", "
1489                    << "None, " << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", "
1490                    << "None, "
1491                    << fitter->second->wmean(
1492                           MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 30.)
1493                    << ", "
1494                    << "None, "
1495                    << fitter->second->wmean(
1496                           MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 20.)
1497                    << ", "
1498                    << "None, "
1499                    << fitter->second->wmean(
1500                           MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 15.)
1501                    << ", "
1502                    << "None, "
1503                    << fitter->second->wmean(
1504                           MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 10.)
1505                    << ", "
1506                    << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", "
1507                    << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", "
1508                    << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", "
1509                    << "None, " << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", "
1510                    << "None)" << std::endl;
1511 
1512             std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
1513             namesimple_x << cname << "_simple_x";
1514             namesimple_dxdz << cname << "_simple_dxdz";
1515             nameweighted_x << cname << "_weighted_x";
1516             nameweighted_dxdz << cname << "_weighted_dxdz";
1517 
1518             fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResid, 10.);
1519             fitter->second->plotsimple(
1520                 namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResSlope, 1000.);
1521 
1522             fitter->second->plotweighted(nameweighted_x.str(),
1523                                          &rootDirectory,
1524                                          MuonResiduals6DOFrphiFitter::kResid,
1525                                          MuonResiduals6DOFrphiFitter::kRedChi2,
1526                                          10.);
1527             fitter->second->plotweighted(nameweighted_dxdz.str(),
1528                                          &rootDirectory,
1529                                          MuonResiduals6DOFrphiFitter::kResSlope,
1530                                          MuonResiduals6DOFrphiFitter::kRedChi2,
1531                                          1000.);
1532           }
1533 
1534           if (successful_fit) {
1535             if (align_x)
1536               params[paramIndex[0]] = deltax_value;
1537             if (align_y)
1538               params[paramIndex[1]] = deltay_value;
1539             if (align_z)
1540               params[paramIndex[2]] = deltaz_value;
1541             if (align_phix)
1542               params[paramIndex[3]] = deltaphix_value;
1543             if (align_phiy)
1544               params[paramIndex[4]] = deltaphiy_value;
1545             if (align_phiz)
1546               params[paramIndex[5]] = deltaphiz_value;
1547           }
1548         }  // end if 6DOFrphi
1549 
1550         if (successful_fit) {
1551           align::Alignables oneortwo;
1552           oneortwo.push_back(ali);
1553           if (thisali != ali)
1554             oneortwo.push_back(thisali);
1555           m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 0., 0.);
1556         } else {
1557           if (m_debug)
1558             std::cout << "MINUIT fit failed!" << std::endl;
1559           if (writeReport) {
1560             report << "reports[-1].status = \"MINUITFAIL\"" << std::endl;
1561           }
1562 
1563           for (int i = 0; i < numParams; i++)
1564             cov[i][i] = 1000.;
1565 
1566           align::Alignables oneortwo;
1567           oneortwo.push_back(ali);
1568           if (thisali != ali)
1569             oneortwo.push_back(thisali);
1570           m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
1571         }
1572       } else {  // too few hits
1573         if (m_debug)
1574           std::cout << "Too few hits!" << std::endl;
1575         if (writeReport) {
1576           report << "reports[-1].status = \"TOOFEWHITS\"" << std::endl;
1577         }
1578 
1579         for (int i = 0; i < numParams; i++)
1580           cov[i][i] = 1000.;
1581 
1582         align::Alignables oneortwo;
1583         oneortwo.push_back(ali);
1584         if (thisali != ali)
1585           oneortwo.push_back(thisali);
1586         m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
1587       }
1588 
1589       AlignmentParameters* parnew = ali->alignmentParameters()->cloneFromSelected(params, cov);
1590       ali->setAlignmentParameters(parnew);
1591       m_alignmentParameterStore->applyParameters(ali);
1592       ali->alignmentParameters()->setValid(true);
1593 
1594       if (m_debug)
1595         std::cout << cname << " fittime= " << stop_watch.CpuTime() << " sec" << std::endl;
1596     }  // end we have a fitter for this alignable
1597 
1598     if (writeReport)
1599       report << std::endl;
1600 
1601   }  // end loop over alignables
1602 
1603   if (writeReport)
1604     report.close();
1605 }
1606 
1607 void MuonAlignmentFromReference::readTmpFiles() {
1608   for (std::vector<std::string>::const_iterator fileName = m_readTemporaryFiles.begin();
1609        fileName != m_readTemporaryFiles.end();
1610        ++fileName) {
1611     FILE* file;
1612     int size;
1613     file = fopen((*fileName).c_str(), "r");
1614     if (file == nullptr)
1615       throw cms::Exception("MuonAlignmentFromReference")
1616           << "file \"" << *fileName << "\" can't be opened (doesn't exist?)" << std::endl;
1617 
1618     fread(&size, sizeof(int), 1, file);
1619     if (int(m_indexes.size()) != size)
1620       throw cms::Exception("MuonAlignmentFromReference")
1621           << "file \"" << *fileName << "\" has " << size << " fitters, but this job has " << m_indexes.size()
1622           << " fitters (probably corresponds to the wrong alignment job)" << std::endl;
1623 
1624     int i = 0;
1625     for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i) {
1626       MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1627       unsigned int index_toread;
1628       fread(&index_toread, sizeof(unsigned int), 1, file);
1629       if (*index != index_toread)
1630         throw cms::Exception("MuonAlignmentFromReference")
1631             << "file \"" << *fileName << "\" has index " << index_toread << " at position " << i
1632             << ", but this job is expecting " << *index << " (probably corresponds to the wrong alignment job)"
1633             << std::endl;
1634       fitter->read(file, i);
1635     }
1636 
1637     fclose(file);
1638   }
1639 }
1640 
1641 void MuonAlignmentFromReference::writeTmpFiles() {
1642   FILE* file;
1643   file = fopen(m_writeTemporaryFile.c_str(), "w");
1644   int size = m_indexes.size();
1645   fwrite(&size, sizeof(int), 1, file);
1646 
1647   int i = 0;
1648   for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i) {
1649     MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1650     unsigned int index_towrite = *index;
1651     fwrite(&index_towrite, sizeof(unsigned int), 1, file);
1652     fitter->write(file, i);
1653   }
1654 
1655   fclose(file);
1656 }
1657 
1658 void MuonAlignmentFromReference::correctBField() {
1659   bool m_debug = false;
1660 
1661   for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1662     if (m_debug)
1663       std::cout << "correcting B in " << chamberPrettyNameFromId(*index) << std::endl;
1664     MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1665     fitter->correctBField();
1666   }
1667 }
1668 
1669 void MuonAlignmentFromReference::fiducialCuts() {
1670   for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1671     if (m_debug)
1672       std::cout << "applying fiducial cuts in " << chamberPrettyNameFromId(*index) << std::endl;
1673     MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1674     fitter->fiducialCuts();
1675   }
1676 }
1677 
1678 void MuonAlignmentFromReference::eraseNotSelectedResiduals() {
1679   for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1680     if (m_debug)
1681       std::cout << "erasing in " << chamberPrettyNameFromId(*index) << std::endl;
1682     MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1683     fitter->eraseNotSelectedResiduals();
1684   }
1685 }
1686 
1687 void MuonAlignmentFromReference::selectResidualsPeaks() {
1688   // should not be called with negative peakNSigma
1689   assert(m_peakNSigma > 0.);
1690 
1691   for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1692     MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1693 
1694     int nvar = 2;
1695     int vars_index[10] = {0, 1};
1696     if (fitter->type() == MuonResidualsFitter::k5DOF) {
1697       if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 ||
1698           fitter->useRes() == MuonResidualsFitter::k1010) {
1699         nvar = 2;
1700         vars_index[0] = MuonResiduals5DOFFitter::kResid;
1701         vars_index[1] = MuonResiduals5DOFFitter::kResSlope;
1702       } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
1703         nvar = 1;
1704         vars_index[0] = MuonResiduals5DOFFitter::kResid;
1705       } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
1706         nvar = 1;
1707         vars_index[0] = MuonResiduals5DOFFitter::kResSlope;
1708       }
1709     } else if (fitter->type() == MuonResidualsFitter::k6DOF) {
1710       if (fitter->useRes() == MuonResidualsFitter::k1111) {
1711         nvar = 4;
1712         vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1713         vars_index[1] = MuonResiduals6DOFFitter::kResidY;
1714         vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
1715         vars_index[3] = MuonResiduals6DOFFitter::kResSlopeY;
1716       } else if (fitter->useRes() == MuonResidualsFitter::k1110) {
1717         nvar = 3;
1718         vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1719         vars_index[1] = MuonResiduals6DOFFitter::kResidY;
1720         vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
1721       } else if (fitter->useRes() == MuonResidualsFitter::k1010) {
1722         nvar = 2;
1723         vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1724         vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
1725       } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
1726         nvar = 2;
1727         vars_index[0] = MuonResiduals6DOFFitter::kResidX;
1728         vars_index[1] = MuonResiduals6DOFFitter::kResidY;
1729       } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
1730         nvar = 1;
1731         vars_index[0] = MuonResiduals6DOFFitter::kResSlopeX;
1732       }
1733     } else if (fitter->type() == MuonResidualsFitter::k6DOFrphi) {
1734       if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 ||
1735           fitter->useRes() == MuonResidualsFitter::k1010) {
1736         nvar = 2;
1737         vars_index[0] = MuonResiduals6DOFrphiFitter::kResid;
1738         vars_index[1] = MuonResiduals6DOFrphiFitter::kResSlope;
1739       } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
1740         nvar = 1;
1741         vars_index[0] = MuonResiduals6DOFrphiFitter::kResid;
1742       } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
1743         nvar = 1;
1744         vars_index[0] = MuonResiduals6DOFrphiFitter::kResSlope;
1745       }
1746     } else
1747       assert(false);
1748 
1749     if (m_debug)
1750       std::cout << "selecting in " << chamberPrettyNameFromId(*index) << std::endl;
1751 
1752     fitter->selectPeakResiduals(m_peakNSigma, nvar, vars_index);
1753   }
1754 }
1755 
1756 std::string MuonAlignmentFromReference::chamberPrettyNameFromId(unsigned int idx) {
1757   DetId id(idx);
1758   char cname[40];
1759   if (id.subdetId() == MuonSubdetId::DT) {
1760     DTChamberId chamberId(id.rawId());
1761     sprintf(cname, "MB%+d/%d/%02d", chamberId.wheel(), chamberId.station(), chamberId.sector());
1762   } else if (id.subdetId() == MuonSubdetId::CSC) {
1763     CSCDetId chamberId(id.rawId());
1764     sprintf(cname,
1765             "ME%s%d/%d/%02d",
1766             (chamberId.endcap() == 1 ? "+" : "-"),
1767             chamberId.station(),
1768             chamberId.ring(),
1769             chamberId.chamber());
1770   }
1771   return std::string(cname);
1772 }
1773 
1774 void MuonAlignmentFromReference::fillNtuple() {
1775   // WARNING: does not support two bin option!!!
1776 
1777   for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index) {
1778     DetId detid(*index);
1779     if (detid.det() != DetId::Muon || !(detid.subdetId() == MuonSubdetId::DT || detid.subdetId() == MuonSubdetId::CSC))
1780       assert(false);
1781 
1782     if (detid.subdetId() == MuonSubdetId::DT) {
1783       m_tree_row.is_dt = (Bool_t) true;
1784       DTChamberId id(*index);
1785       m_tree_row.is_plus = (Bool_t) true;
1786       m_tree_row.station = (UChar_t)id.station();
1787       m_tree_row.ring_wheel = (Char_t)id.wheel();
1788       m_tree_row.sector = (UChar_t)id.sector();
1789     } else {
1790       m_tree_row.is_dt = (Bool_t) false;
1791       CSCDetId id(*index);
1792       m_tree_row.is_plus = (Bool_t)(id.endcap() == 1);
1793       m_tree_row.station = (UChar_t)id.station();
1794       m_tree_row.ring_wheel = (Char_t)id.ring();
1795       m_tree_row.sector = (UChar_t)id.chamber();
1796     }
1797 
1798     MuonResidualsTwoBin* fitter = m_fitterOrder[*index];
1799 
1800     std::vector<double*>::const_iterator residual = fitter->residualsPos_begin();
1801     std::vector<bool>::const_iterator residual_ok = fitter->residualsPos_ok_begin();
1802     for (; residual != fitter->residualsPos_end(); ++residual, ++residual_ok) {
1803       if (fitter->type() == MuonResidualsFitter::k5DOF || fitter->type() == MuonResidualsFitter::k6DOFrphi) {
1804         m_tree_row.res_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kResid];
1805         m_tree_row.res_y = (Float_t)0.;
1806         m_tree_row.res_slope_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kResSlope];
1807         m_tree_row.res_slope_y = (Float_t)0.;
1808         m_tree_row.pos_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPositionX];
1809         m_tree_row.pos_y = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPositionY];
1810         m_tree_row.angle_x = (Float_t)(*residual)[MuonResiduals5DOFFitter::kAngleX];
1811         m_tree_row.angle_y = (Float_t)(*residual)[MuonResiduals5DOFFitter::kAngleY];
1812         m_tree_row.pz = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPz];
1813         m_tree_row.pt = (Float_t)(*residual)[MuonResiduals5DOFFitter::kPt];
1814         m_tree_row.q = (Char_t)(*residual)[MuonResiduals5DOFFitter::kCharge];
1815         m_tree_row.select = (Bool_t)*residual_ok;
1816       } else if (fitter->type() == MuonResidualsFitter::k6DOF) {
1817         m_tree_row.res_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResidX];
1818         m_tree_row.res_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResidY];
1819         m_tree_row.res_slope_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResSlopeX];
1820         m_tree_row.res_slope_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kResSlopeY];
1821         m_tree_row.pos_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPositionX];
1822         m_tree_row.pos_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPositionY];
1823         m_tree_row.angle_x = (Float_t)(*residual)[MuonResiduals6DOFFitter::kAngleX];
1824         m_tree_row.angle_y = (Float_t)(*residual)[MuonResiduals6DOFFitter::kAngleY];
1825         m_tree_row.pz = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPz];
1826         m_tree_row.pt = (Float_t)(*residual)[MuonResiduals6DOFFitter::kPt];
1827         m_tree_row.q = (Char_t)(*residual)[MuonResiduals6DOFFitter::kCharge];
1828         m_tree_row.select = (Bool_t)*residual_ok;
1829       } else
1830         assert(false);
1831 
1832       m_ttree->Fill();
1833     }
1834   }
1835 }
1836 
1837 void MuonAlignmentFromReference::parseReference(align::Alignables& reference,
1838                                                 const align::Alignables& all_DT_chambers,
1839                                                 const align::Alignables& all_CSC_chambers) {
1840   std::map<Alignable*, bool> already_seen;
1841 
1842   for (std::vector<std::string>::const_iterator name = m_reference.begin(); name != m_reference.end(); ++name) {
1843     bool parsing_error = false;
1844 
1845     bool barrel = (name->substr(0, 2) == std::string("MB"));
1846     bool endcap = (name->substr(0, 2) == std::string("ME"));
1847     if (!barrel && !endcap)
1848       parsing_error = true;
1849 
1850     if (!parsing_error && barrel) {
1851       int index = 2;
1852       if (name->substr(index, 1) == std::string(" "))
1853         index++;
1854 
1855       bool plus = true;
1856       if (name->substr(index, 1) == std::string("+")) {
1857         plus = true;
1858         index++;
1859       } else if (name->substr(index, 1) == std::string("-")) {
1860         plus = false;
1861         index++;
1862       } else if (numeric(name->substr(index, 1))) {
1863       } else
1864         parsing_error = true;
1865 
1866       int wheel = 0;
1867       bool wheel_digit = false;
1868       while (!parsing_error && numeric(name->substr(index, 1))) {
1869         wheel *= 10;
1870         wheel += number(name->substr(index, 1));
1871         wheel_digit = true;
1872         index++;
1873       }
1874       if (!plus)
1875         wheel *= -1;
1876       if (!wheel_digit)
1877         parsing_error = true;
1878 
1879       if (name->substr(index, 1) != std::string(" "))
1880         parsing_error = true;
1881       index++;
1882 
1883       int station = 0;
1884       bool station_digit = false;
1885       while (!parsing_error && numeric(name->substr(index, 1))) {
1886         station *= 10;
1887         station += number(name->substr(index, 1));
1888         station_digit = true;
1889         index++;
1890       }
1891       if (!station_digit)
1892         parsing_error = true;
1893 
1894       if (name->substr(index, 1) != std::string(" "))
1895         parsing_error = true;
1896       index++;
1897 
1898       int sector = 0;
1899       bool sector_digit = false;
1900       while (!parsing_error && numeric(name->substr(index, 1))) {
1901         sector *= 10;
1902         sector += number(name->substr(index, 1));
1903         sector_digit = true;
1904         index++;
1905       }
1906       if (!sector_digit)
1907         parsing_error = true;
1908 
1909       if (!parsing_error) {
1910         bool no_such_chamber = false;
1911 
1912         if (wheel < -2 || wheel > 2)
1913           no_such_chamber = true;
1914         if (station < 1 || station > 4)
1915           no_such_chamber = true;
1916         if (station == 4 && (sector < 1 || sector > 14))
1917           no_such_chamber = true;
1918         if (station < 4 && (sector < 1 || sector > 12))
1919           no_such_chamber = true;
1920 
1921         if (no_such_chamber)
1922           throw cms::Exception("MuonAlignmentFromReference")
1923               << "reference chamber doesn't exist: " << (*name) << std::endl;
1924 
1925         DTChamberId id(wheel, station, sector);
1926         for (const auto& ali : all_DT_chambers) {
1927           if (ali->geomDetId().rawId() == id.rawId()) {
1928             std::map<Alignable*, bool>::const_iterator trial = already_seen.find(ali);
1929             if (trial == already_seen.end()) {
1930               reference.push_back(ali);
1931               already_seen[ali] = true;
1932             }
1933           }
1934         }
1935       }  // if (!parsing_error)
1936     }
1937 
1938     if (!parsing_error && endcap) {
1939       int index = 2;
1940       if (name->substr(index, 1) == std::string(" "))
1941         index++;
1942 
1943       bool plus = true;
1944       if (name->substr(index, 1) == std::string("+")) {
1945         plus = true;
1946         index++;
1947       } else if (name->substr(index, 1) == std::string("-")) {
1948         plus = false;
1949         index++;
1950       } else if (numeric(name->substr(index, 1))) {
1951       } else
1952         parsing_error = true;
1953 
1954       int station = 0;
1955       bool station_digit = false;
1956       while (!parsing_error && numeric(name->substr(index, 1))) {
1957         station *= 10;
1958         station += number(name->substr(index, 1));
1959         station_digit = true;
1960         index++;
1961       }
1962       if (!plus)
1963         station *= -1;
1964       if (!station_digit)
1965         parsing_error = true;
1966 
1967       if (name->substr(index, 1) != std::string("/"))
1968         parsing_error = true;
1969       index++;
1970 
1971       int ring = 0;
1972       bool ring_digit = false;
1973       while (!parsing_error && numeric(name->substr(index, 1))) {
1974         ring *= 10;
1975         ring += number(name->substr(index, 1));
1976         ring_digit = true;
1977         index++;
1978       }
1979       if (!ring_digit)
1980         parsing_error = true;
1981 
1982       if (name->substr(index, 1) != std::string(" "))
1983         parsing_error = true;
1984       index++;
1985 
1986       int chamber = 0;
1987       bool chamber_digit = false;
1988       while (!parsing_error && numeric(name->substr(index, 1))) {
1989         chamber *= 10;
1990         chamber += number(name->substr(index, 1));
1991         chamber_digit = true;
1992         index++;
1993       }
1994       if (!chamber_digit)
1995         parsing_error = true;
1996 
1997       if (!parsing_error) {
1998         bool no_such_chamber = false;
1999 
2000         int endcap = (station > 0 ? 1 : 2);
2001         station = abs(station);
2002         if (station < 1 || station > 4)
2003           no_such_chamber = true;
2004         if (station == 1 && (ring < 1 || ring > 4))
2005           no_such_chamber = true;
2006         if (station > 1 && (ring < 1 || ring > 2))
2007           no_such_chamber = true;
2008         if (station == 1 && (chamber < 1 || chamber > 36))
2009           no_such_chamber = true;
2010         if (station > 1 && ring == 1 && (chamber < 1 || chamber > 18))
2011           no_such_chamber = true;
2012         if (station > 1 && ring == 2 && (chamber < 1 || chamber > 36))
2013           no_such_chamber = true;
2014 
2015         if (no_such_chamber)
2016           throw cms::Exception("MuonAlignmentFromReference")
2017               << "reference chamber doesn't exist: " << (*name) << std::endl;
2018 
2019         CSCDetId id(endcap, station, ring, chamber);
2020         for (const auto& ali : all_CSC_chambers) {
2021           if (ali->geomDetId().rawId() == id.rawId()) {
2022             std::map<Alignable*, bool>::const_iterator trial = already_seen.find(ali);
2023             if (trial == already_seen.end()) {
2024               reference.push_back(ali);
2025               already_seen[ali] = true;
2026             }
2027           }
2028         }
2029       }  // if (!parsing_error)
2030     }  // endcap
2031 
2032     if (parsing_error)
2033       throw cms::Exception("MuonAlignmentFromReference")
2034           << "reference chamber name is malformed: " << (*name) << std::endl;
2035   }
2036 }
2037 
2038 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
2039 DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory, MuonAlignmentFromReference, "MuonAlignmentFromReference");
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!