Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-20 10:41:40

0001 #include "Alignment/MuonAlignmentAlgorithms/plugins/CSCOverlapsAlignmentAlgorithm.h"
0002 
0003 CSCOverlapsAlignmentAlgorithm::CSCOverlapsAlignmentAlgorithm(const edm::ParameterSet& iConfig,
0004                                                              edm::ConsumesCollector& iC)
0005     : AlignmentAlgorithmBase(iConfig, iC),
0006       m_minHitsPerChamber(iConfig.getParameter<int>("minHitsPerChamber")),
0007       m_maxdrdz(iConfig.getParameter<double>("maxdrdz")),
0008       m_fiducial(iConfig.getParameter<bool>("fiducial")),
0009       m_useHitWeights(iConfig.getParameter<bool>("useHitWeights")),
0010       m_slopeFromTrackRefit(iConfig.getParameter<bool>("slopeFromTrackRefit")),
0011       m_minStationsInTrackRefits(iConfig.getParameter<int>("minStationsInTrackRefits")),
0012       m_truncateSlopeResid(iConfig.getParameter<double>("truncateSlopeResid")),
0013       m_truncateOffsetResid(iConfig.getParameter<double>("truncateOffsetResid")),
0014       m_combineME11(iConfig.getParameter<bool>("combineME11")),
0015       m_useTrackWeights(iConfig.getParameter<bool>("useTrackWeights")),
0016       m_errorFromRMS(iConfig.getParameter<bool>("errorFromRMS")),
0017       m_minTracksPerOverlap(iConfig.getParameter<int>("minTracksPerOverlap")),
0018       m_makeHistograms(iConfig.getParameter<bool>("makeHistograms")),
0019       m_cscGeometryToken(iC.esConsumes<edm::Transition::BeginRun>()),
0020       m_propToken(iC.esConsumes(edm::ESInputTag(
0021           "",
0022           m_slopeFromTrackRefit
0023               ? iConfig.getParameter<edm::ParameterSet>("TrackTransformer").getParameter<std::string>("Propagator")
0024               : std::string("")))),
0025       m_tthbToken(iC.esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0026       m_mode_string(iConfig.getParameter<std::string>("mode")),
0027       m_reportFileName(iConfig.getParameter<std::string>("reportFileName")),
0028       m_minP(iConfig.getParameter<double>("minP")),
0029       m_maxRedChi2(iConfig.getParameter<double>("maxRedChi2")),
0030       m_writeTemporaryFile(iConfig.getParameter<std::string>("writeTemporaryFile")),
0031       m_readTemporaryFiles(iConfig.getParameter<std::vector<std::string> >("readTemporaryFiles")),
0032       m_doAlignment(iConfig.getParameter<bool>("doAlignment")) {
0033   if (m_mode_string == std::string("phiy"))
0034     m_mode = CSCPairResidualsConstraint::kModePhiy;
0035   else if (m_mode_string == std::string("phipos"))
0036     m_mode = CSCPairResidualsConstraint::kModePhiPos;
0037   else if (m_mode_string == std::string("phiz"))
0038     m_mode = CSCPairResidualsConstraint::kModePhiz;
0039   else if (m_mode_string == std::string("radius"))
0040     m_mode = CSCPairResidualsConstraint::kModeRadius;
0041   else
0042     throw cms::Exception("BadConfig") << "mode must be one of \"phiy\", \"phipos\", \"phiz\", \"radius\"" << std::endl;
0043 
0044   std::vector<edm::ParameterSet> fitters = iConfig.getParameter<std::vector<edm::ParameterSet> >("fitters");
0045   for (std::vector<edm::ParameterSet>::const_iterator fitter = fitters.begin(); fitter != fitters.end(); ++fitter) {
0046     m_fitters.push_back(CSCChamberFitter(*fitter, m_residualsConstraints));
0047   }
0048 
0049   for (std::vector<CSCPairResidualsConstraint*>::const_iterator residualsConstraint = m_residualsConstraints.begin();
0050        residualsConstraint != m_residualsConstraints.end();
0051        ++residualsConstraint) {
0052     (*residualsConstraint)->configure(this);
0053     m_quickChamberLookup[std::pair<CSCDetId, CSCDetId>((*residualsConstraint)->id_i(), (*residualsConstraint)->id_j())] =
0054         *residualsConstraint;
0055   }
0056 
0057   if (m_slopeFromTrackRefit) {
0058     m_trackTransformer = new TrackTransformer(iConfig.getParameter<edm::ParameterSet>("TrackTransformer"), iC);
0059   } else {
0060     m_trackTransformer = nullptr;
0061   }
0062 
0063   m_propagatorPointer = nullptr;
0064 
0065   if (m_makeHistograms) {
0066     edm::Service<TFileService> tFileService;
0067     m_histP10 = tFileService->make<TH1F>("P10", "", 100, 0, 10);
0068     m_histP100 = tFileService->make<TH1F>("P100", "", 100, 0, 100);
0069     m_histP1000 = tFileService->make<TH1F>("P1000", "", 100, 0, 1000);
0070 
0071     m_hitsPerChamber = tFileService->make<TH1F>("hitsPerChamber", "", 10, -0.5, 9.5);
0072 
0073     m_fiducial_ME11 = tFileService->make<TProfile>("fiducial_ME11", "", 100, 0.075, 0.100);
0074     m_fiducial_ME12 = tFileService->make<TProfile>("fiducial_ME12", "", 100, 0.080, 0.105);
0075     m_fiducial_MEx1 = tFileService->make<TProfile>("fiducial_MEx1", "", 100, 0.160, 0.210);
0076     m_fiducial_MEx2 = tFileService->make<TProfile>("fiducial_MEx2", "", 100, 0.080, 0.105);
0077 
0078     m_slope = tFileService->make<TH1F>("slope", "", 100, -0.5, 0.5);
0079     m_slope_MEp4 = tFileService->make<TH1F>("slope_MEp4", "", 100, -0.5, 0.5);
0080     m_slope_MEp3 = tFileService->make<TH1F>("slope_MEp3", "", 100, -0.5, 0.5);
0081     m_slope_MEp2 = tFileService->make<TH1F>("slope_MEp2", "", 100, -0.5, 0.5);
0082     m_slope_MEp1 = tFileService->make<TH1F>("slope_MEp1", "", 100, -0.5, 0.5);
0083     m_slope_MEm1 = tFileService->make<TH1F>("slope_MEm1", "", 100, -0.5, 0.5);
0084     m_slope_MEm2 = tFileService->make<TH1F>("slope_MEm2", "", 100, -0.5, 0.5);
0085     m_slope_MEm3 = tFileService->make<TH1F>("slope_MEm3", "", 100, -0.5, 0.5);
0086     m_slope_MEm4 = tFileService->make<TH1F>("slope_MEm4", "", 100, -0.5, 0.5);
0087 
0088     m_slopeResiduals = tFileService->make<TH1F>("slopeResiduals", "mrad", 300, -30., 30.);
0089     m_slopeResiduals_weighted = tFileService->make<TH1F>("slopeResiduals_weighted", "mrad", 300, -30., 30.);
0090     m_slopeResiduals_normalized = tFileService->make<TH1F>("slopeResiduals_normalized", "", 200, -20., 20.);
0091     m_offsetResiduals = tFileService->make<TH1F>("offsetResiduals", "mm", 300, -30., 30.);
0092     m_offsetResiduals_weighted = tFileService->make<TH1F>("offsetResiduals_weighted", "mm", 300, -30., 30.);
0093     m_offsetResiduals_normalized = tFileService->make<TH1F>("offsetResiduals_normalized", "", 200, -20., 20.);
0094 
0095     m_drdz = tFileService->make<TH1F>("drdz", "", 100, -0.5, 0.5);
0096 
0097     m_occupancy = tFileService->make<TH2F>("occupancy", "", 36, 1, 37, 20, 1, 21);
0098     for (int i = 1; i <= 36; i++) {
0099       std::stringstream pairname;
0100       pairname << i << "-";
0101       if (i + 1 == 37)
0102         pairname << 1;
0103       else
0104         pairname << (i + 1);
0105       m_occupancy->GetXaxis()->SetBinLabel(i, pairname.str().c_str());
0106     }
0107     m_occupancy->GetYaxis()->SetBinLabel(1, "ME-4/2");
0108     m_occupancy->GetYaxis()->SetBinLabel(2, "ME-4/1");
0109     m_occupancy->GetYaxis()->SetBinLabel(3, "ME-3/2");
0110     m_occupancy->GetYaxis()->SetBinLabel(4, "ME-3/1");
0111     m_occupancy->GetYaxis()->SetBinLabel(5, "ME-2/2");
0112     m_occupancy->GetYaxis()->SetBinLabel(6, "ME-2/1");
0113     m_occupancy->GetYaxis()->SetBinLabel(7, "ME-1/3");
0114     m_occupancy->GetYaxis()->SetBinLabel(8, "ME-1/2");
0115     if (!m_combineME11) {
0116       m_occupancy->GetYaxis()->SetBinLabel(9, "ME-1/1b");
0117       m_occupancy->GetYaxis()->SetBinLabel(10, "ME-1/1a");
0118       m_occupancy->GetYaxis()->SetBinLabel(11, "ME+1/1a");
0119       m_occupancy->GetYaxis()->SetBinLabel(12, "ME+1/1b");
0120     } else {
0121       m_occupancy->GetYaxis()->SetBinLabel(9, "ME-1/1");
0122       m_occupancy->GetYaxis()->SetBinLabel(10, "");
0123       m_occupancy->GetYaxis()->SetBinLabel(11, "");
0124       m_occupancy->GetYaxis()->SetBinLabel(12, "ME+1/1");
0125     }
0126     m_occupancy->GetYaxis()->SetBinLabel(13, "ME+1/2");
0127     m_occupancy->GetYaxis()->SetBinLabel(14, "ME+1/3");
0128     m_occupancy->GetYaxis()->SetBinLabel(15, "ME+2/1");
0129     m_occupancy->GetYaxis()->SetBinLabel(16, "ME+2/2");
0130     m_occupancy->GetYaxis()->SetBinLabel(17, "ME+3/1");
0131     m_occupancy->GetYaxis()->SetBinLabel(18, "ME+3/2");
0132     m_occupancy->GetYaxis()->SetBinLabel(19, "ME+4/1");
0133     m_occupancy->GetYaxis()->SetBinLabel(20, "ME+4/2");
0134 
0135     m_XYpos_mep1 = tFileService->make<TH2F>("XYpos_mep1", "Positions: ME+1", 140, -700., 700., 140, -700., 700.);
0136     m_XYpos_mep2 = tFileService->make<TH2F>("XYpos_mep2", "Positions: ME+2", 140, -700., 700., 140, -700., 700.);
0137     m_XYpos_mep3 = tFileService->make<TH2F>("XYpos_mep3", "Positions: ME+3", 140, -700., 700., 140, -700., 700.);
0138     m_XYpos_mep4 = tFileService->make<TH2F>("XYpos_mep4", "Positions: ME+4", 140, -700., 700., 140, -700., 700.);
0139     m_XYpos_mem1 = tFileService->make<TH2F>("XYpos_mem1", "Positions: ME-1", 140, -700., 700., 140, -700., 700.);
0140     m_XYpos_mem2 = tFileService->make<TH2F>("XYpos_mem2", "Positions: ME-2", 140, -700., 700., 140, -700., 700.);
0141     m_XYpos_mem3 = tFileService->make<TH2F>("XYpos_mem3", "Positions: ME-3", 140, -700., 700., 140, -700., 700.);
0142     m_XYpos_mem4 = tFileService->make<TH2F>("XYpos_mem4", "Positions: ME-4", 140, -700., 700., 140, -700., 700.);
0143     m_RPhipos_mep1 = tFileService->make<TH2F>("RPhipos_mep1", "Positions: ME+1", 144, -M_PI, M_PI, 21, 0., 700.);
0144     m_RPhipos_mep2 = tFileService->make<TH2F>("RPhipos_mep2", "Positions: ME+2", 144, -M_PI, M_PI, 21, 0., 700.);
0145     m_RPhipos_mep3 = tFileService->make<TH2F>("RPhipos_mep3", "Positions: ME+3", 144, -M_PI, M_PI, 21, 0., 700.);
0146     m_RPhipos_mep4 = tFileService->make<TH2F>("RPhipos_mep4", "Positions: ME+4", 144, -M_PI, M_PI, 21, 0., 700.);
0147     m_RPhipos_mem1 = tFileService->make<TH2F>("RPhipos_mem1", "Positions: ME-1", 144, -M_PI, M_PI, 21, 0., 700.);
0148     m_RPhipos_mem2 = tFileService->make<TH2F>("RPhipos_mem2", "Positions: ME-2", 144, -M_PI, M_PI, 21, 0., 700.);
0149     m_RPhipos_mem3 = tFileService->make<TH2F>("RPhipos_mem3", "Positions: ME-3", 144, -M_PI, M_PI, 21, 0., 700.);
0150     m_RPhipos_mem4 = tFileService->make<TH2F>("RPhipos_mem4", "Positions: ME-4", 144, -M_PI, M_PI, 21, 0., 700.);
0151   } else {
0152     m_histP10 = nullptr;
0153     m_histP100 = nullptr;
0154     m_histP1000 = nullptr;
0155     m_hitsPerChamber = nullptr;
0156     m_fiducial_ME11 = nullptr;
0157     m_fiducial_ME12 = nullptr;
0158     m_fiducial_MEx1 = nullptr;
0159     m_fiducial_MEx2 = nullptr;
0160     m_slope = nullptr;
0161     m_slope_MEp4 = nullptr;
0162     m_slope_MEp3 = nullptr;
0163     m_slope_MEp2 = nullptr;
0164     m_slope_MEp1 = nullptr;
0165     m_slope_MEm1 = nullptr;
0166     m_slope_MEm2 = nullptr;
0167     m_slope_MEm3 = nullptr;
0168     m_slope_MEm4 = nullptr;
0169     m_slopeResiduals = nullptr;
0170     m_slopeResiduals_weighted = nullptr;
0171     m_slopeResiduals_normalized = nullptr;
0172     m_offsetResiduals = nullptr;
0173     m_offsetResiduals_weighted = nullptr;
0174     m_offsetResiduals_normalized = nullptr;
0175     m_drdz = nullptr;
0176     m_occupancy = nullptr;
0177     m_XYpos_mep1 = nullptr;
0178     m_XYpos_mep2 = nullptr;
0179     m_XYpos_mep3 = nullptr;
0180     m_XYpos_mep4 = nullptr;
0181     m_XYpos_mem1 = nullptr;
0182     m_XYpos_mem2 = nullptr;
0183     m_XYpos_mem3 = nullptr;
0184     m_XYpos_mem4 = nullptr;
0185     m_RPhipos_mep1 = nullptr;
0186     m_RPhipos_mep2 = nullptr;
0187     m_RPhipos_mep3 = nullptr;
0188     m_RPhipos_mep4 = nullptr;
0189     m_RPhipos_mem1 = nullptr;
0190     m_RPhipos_mem2 = nullptr;
0191     m_RPhipos_mem3 = nullptr;
0192     m_RPhipos_mem4 = nullptr;
0193   }
0194 }
0195 
0196 CSCOverlapsAlignmentAlgorithm::~CSCOverlapsAlignmentAlgorithm() {}
0197 
0198 void CSCOverlapsAlignmentAlgorithm::initialize(const edm::EventSetup& iSetup,
0199                                                AlignableTracker* alignableTracker,
0200                                                AlignableMuon* alignableMuon,
0201                                                AlignableExtras* alignableExtras,
0202                                                AlignmentParameterStore* alignmentParameterStore) {
0203   m_alignmentParameterStore = alignmentParameterStore;
0204   m_alignables = m_alignmentParameterStore->alignables();
0205 
0206   if (alignableTracker == nullptr)
0207     m_alignableNavigator = new AlignableNavigator(alignableMuon);
0208   else
0209     m_alignableNavigator = new AlignableNavigator(alignableTracker, alignableMuon);
0210 
0211   for (const auto& alignable : m_alignables) {
0212     DetId id = alignable->geomDetId();
0213     if (id.det() != DetId::Muon || id.subdetId() != MuonSubdetId::CSC || CSCDetId(id.rawId()).layer() != 0) {
0214       throw cms::Exception("BadConfig") << "Only CSC chambers may be alignable" << std::endl;
0215     }
0216 
0217     std::vector<bool> selector = alignable->alignmentParameters()->selector();
0218     for (std::vector<bool>::const_iterator i = selector.begin(); i != selector.end(); ++i) {
0219       if (!(*i))
0220         throw cms::Exception("BadConfig") << "All selector strings should be \"111111\"" << std::endl;
0221     }
0222   }
0223 
0224   const CSCGeometry* cscGeometry = &iSetup.getData(m_cscGeometryToken);
0225 
0226   for (std::vector<CSCPairResidualsConstraint*>::const_iterator residualsConstraint = m_residualsConstraints.begin();
0227        residualsConstraint != m_residualsConstraints.end();
0228        ++residualsConstraint) {
0229     (*residualsConstraint)->setZplane(cscGeometry);
0230   }
0231 
0232   if (!m_readTemporaryFiles.empty()) {
0233     std::vector<std::ifstream*> input;
0234     for (std::vector<std::string>::const_iterator fileName = m_readTemporaryFiles.begin();
0235          fileName != m_readTemporaryFiles.end();
0236          ++fileName) {
0237       input.push_back(new std::ifstream(fileName->c_str()));
0238     }
0239 
0240     for (std::vector<CSCPairResidualsConstraint*>::const_iterator residualsConstraint = m_residualsConstraints.begin();
0241          residualsConstraint != m_residualsConstraints.end();
0242          ++residualsConstraint) {
0243       (*residualsConstraint)->read(input, m_readTemporaryFiles);
0244     }
0245 
0246     for (std::vector<std::ifstream*>::const_iterator file = input.begin(); file != input.end(); ++file) {
0247       delete (*file);
0248     }
0249   }
0250 }
0251 
0252 void CSCOverlapsAlignmentAlgorithm::run(const edm::EventSetup& iSetup, const EventInfo& eventInfo) {
0253   edm::ESHandle<Propagator> propagator;
0254   if (m_slopeFromTrackRefit) {
0255     iSetup.getHandle(m_propToken);
0256     if (m_propagatorPointer != &*propagator) {
0257       m_propagatorPointer = &*propagator;
0258 
0259       for (std::vector<CSCPairResidualsConstraint*>::const_iterator residualsConstraint =
0260                m_residualsConstraints.begin();
0261            residualsConstraint != m_residualsConstraints.end();
0262            ++residualsConstraint) {
0263         (*residualsConstraint)->setPropagator(m_propagatorPointer);
0264       }
0265     }
0266   }
0267 
0268   const TransientTrackBuilder* transientTrackBuilder = &iSetup.getData(m_tthbToken);
0269 
0270   if (m_trackTransformer != nullptr)
0271     m_trackTransformer->setServices(iSetup);
0272 
0273   const ConstTrajTrackPairCollection& trajtracks = eventInfo.trajTrackPairs();
0274   for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
0275        ++trajtrack) {
0276     const Trajectory* traj = (*trajtrack).first;
0277     const reco::Track* track = (*trajtrack).second;
0278 
0279     if (m_makeHistograms) {
0280       m_histP10->Fill(track->p());
0281       m_histP100->Fill(track->p());
0282       m_histP1000->Fill(track->p());
0283     }
0284     if (track->p() >= m_minP) {
0285       std::vector<TrajectoryMeasurement> measurements = traj->measurements();
0286       reco::TransientTrack transientTrack = transientTrackBuilder->build(track);
0287 
0288       std::map<int, std::map<CSCDetId, bool> > stationsToChambers;
0289       for (std::vector<TrajectoryMeasurement>::const_iterator measurement = measurements.begin();
0290            measurement != measurements.end();
0291            ++measurement) {
0292         DetId id = measurement->recHit()->geographicalId();
0293         if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
0294           CSCDetId cscid(id.rawId());
0295           CSCDetId chamberId(cscid.endcap(), cscid.station(), cscid.ring(), cscid.chamber(), 0);
0296           if (m_combineME11 && cscid.station() == 1 && cscid.ring() == 4)
0297             chamberId = CSCDetId(cscid.endcap(), 1, 1, cscid.chamber(), 0);
0298           int station = (cscid.endcap() == 1 ? 1 : -1) * cscid.station();
0299 
0300           if (stationsToChambers.find(station) == stationsToChambers.end())
0301             stationsToChambers[station] = std::map<CSCDetId, bool>();
0302           stationsToChambers[station][chamberId] = true;
0303 
0304           if (m_makeHistograms) {
0305             GlobalPoint pos = measurement->recHit()->globalPosition();
0306             if (cscid.endcap() == 1 && cscid.station() == 1) {
0307               m_XYpos_mep1->Fill(pos.x(), pos.y());
0308               m_RPhipos_mep1->Fill(pos.phi(), pos.perp());
0309             }
0310             if (cscid.endcap() == 1 && cscid.station() == 2) {
0311               m_XYpos_mep2->Fill(pos.x(), pos.y());
0312               m_RPhipos_mep2->Fill(pos.phi(), pos.perp());
0313             }
0314             if (cscid.endcap() == 1 && cscid.station() == 3) {
0315               m_XYpos_mep3->Fill(pos.x(), pos.y());
0316               m_RPhipos_mep3->Fill(pos.phi(), pos.perp());
0317             }
0318             if (cscid.endcap() == 1 && cscid.station() == 4) {
0319               m_XYpos_mep4->Fill(pos.x(), pos.y());
0320               m_RPhipos_mep4->Fill(pos.phi(), pos.perp());
0321             }
0322             if (cscid.endcap() == 2 && cscid.station() == 1) {
0323               m_XYpos_mem1->Fill(pos.x(), pos.y());
0324               m_RPhipos_mem1->Fill(pos.phi(), pos.perp());
0325             }
0326             if (cscid.endcap() == 2 && cscid.station() == 2) {
0327               m_XYpos_mem2->Fill(pos.x(), pos.y());
0328               m_RPhipos_mem2->Fill(pos.phi(), pos.perp());
0329             }
0330             if (cscid.endcap() == 2 && cscid.station() == 3) {
0331               m_XYpos_mem3->Fill(pos.x(), pos.y());
0332               m_RPhipos_mem3->Fill(pos.phi(), pos.perp());
0333             }
0334             if (cscid.endcap() == 2 && cscid.station() == 4) {
0335               m_XYpos_mem4->Fill(pos.x(), pos.y());
0336               m_RPhipos_mem4->Fill(pos.phi(), pos.perp());
0337             }
0338           }
0339         }
0340       }
0341 
0342       std::map<CSCPairResidualsConstraint*, bool> residualsConstraints;
0343       for (std::map<int, std::map<CSCDetId, bool> >::const_iterator iter = stationsToChambers.begin();
0344            iter != stationsToChambers.end();
0345            ++iter) {
0346         for (std::map<CSCDetId, bool>::const_iterator one = iter->second.begin(); one != iter->second.end(); ++one) {
0347           for (std::map<CSCDetId, bool>::const_iterator two = one; two != iter->second.end(); ++two) {
0348             if (one != two) {
0349               std::map<std::pair<CSCDetId, CSCDetId>, CSCPairResidualsConstraint*>::const_iterator quick;
0350 
0351               quick = m_quickChamberLookup.find(std::pair<CSCDetId, CSCDetId>(one->first, two->first));
0352               if (quick != m_quickChamberLookup.end())
0353                 residualsConstraints[quick->second] = true;
0354 
0355               quick = m_quickChamberLookup.find(std::pair<CSCDetId, CSCDetId>(two->first, one->first));
0356               if (quick != m_quickChamberLookup.end())
0357                 residualsConstraints[quick->second] = true;
0358             }
0359           }
0360         }
0361       }
0362 
0363       for (std::map<CSCPairResidualsConstraint*, bool>::const_iterator residualsConstraint =
0364                residualsConstraints.begin();
0365            residualsConstraint != residualsConstraints.end();
0366            ++residualsConstraint) {
0367         residualsConstraint->first->addTrack(measurements, transientTrack, m_trackTransformer);
0368       }
0369     }
0370   }
0371 }
0372 
0373 void CSCOverlapsAlignmentAlgorithm::terminate(const edm::EventSetup& iSetup) {
0374   // write residuals partial fits to temporary files for collection
0375   if (m_writeTemporaryFile != std::string("")) {
0376     std::ofstream output(m_writeTemporaryFile.c_str());
0377     for (std::vector<CSCPairResidualsConstraint*>::const_iterator residualsConstraint = m_residualsConstraints.begin();
0378          residualsConstraint != m_residualsConstraints.end();
0379          ++residualsConstraint) {
0380       (*residualsConstraint)->write(output);
0381     }
0382   }
0383 
0384   // write report for alignment results
0385   if (m_doAlignment) {
0386     std::ofstream report;
0387     bool writeReport = (m_reportFileName != std::string(""));
0388     if (writeReport) {
0389       report.open(m_reportFileName.c_str());
0390       report << "cscReports = []" << std::endl
0391              << std::endl
0392              << "class CSCChamberCorrection:" << std::endl
0393              << "    def __init__(self, name, detid, value):" << std::endl
0394              << "        self.name, self.detid, self.value = name, detid, value" << std::endl
0395              << std::endl
0396              << "class CSCErrorMode:" << std::endl
0397              << "    def __init__(self, error):" << std::endl
0398              << "        self.error = error" << std::endl
0399              << "        self.terms = {}" << std::endl
0400              << "        self.detids = {}" << std::endl
0401              << "    def addTerm(self, name, detid, coefficient):" << std::endl
0402              << "        self.terms[name] = coefficient" << std::endl
0403              << "        self.detids[name] = detid" << std::endl
0404              << std::endl
0405              << "class CSCConstraintResidual:" << std::endl
0406              << "    def __init__(self, i, j, before, uncert, residual, pull):" << std::endl
0407              << "        self.i, self.j, self.before, self.error, self.residual, self.pull = i, j, before, uncert, "
0408                 "residual, pull"
0409              << std::endl
0410              << std::endl
0411              << "class CSCFitterReport:" << std::endl
0412              << "    def __init__(self, name, oldchi2, newchi2):" << std::endl
0413              << "        self.name, self.oldchi2, self.newchi2 = name, oldchi2, newchi2" << std::endl
0414              << "        self.chamberCorrections = []" << std::endl
0415              << "        self.errorModes = []" << std::endl
0416              << "        self.constraintResiduals = []" << std::endl
0417              << std::endl
0418              << "    def addChamberCorrection(self, name, detid, value):" << std::endl
0419              << "        self.chamberCorrections.append(CSCChamberCorrection(name, detid, value))" << std::endl
0420              << std::endl
0421              << "    def addErrorMode(self, error):" << std::endl
0422              << "        self.errorModes.append(CSCErrorMode(error))" << std::endl
0423              << std::endl
0424              << "    def addErrorModeTerm(self, name, detid, coefficient):" << std::endl
0425              << "        self.errorModes[-1].addTerm(name, detid, coefficient)" << std::endl
0426              << std::endl
0427              << "    def addCSCConstraintResidual(self, i, j, before, uncert, residual, pull):" << std::endl
0428              << "        self.constraintResiduals.append(CSCConstraintResidual(i, j, before, uncert, residual, pull))"
0429              << std::endl
0430              << std::endl
0431              << "import re" << std::endl
0432              << "def nameToKey(name):" << std::endl
0433              << "    match = re.match(\"ME([\\+\\-])([1-4])/([1-4])/([0-9]{2})\", name)" << std::endl
0434              << "    if match is None: return None" << std::endl
0435              << "    endcap, station, ring, chamber = match.groups()" << std::endl
0436              << "    if endcap == \"+\": endcap = 1" << std::endl
0437              << "    else: endcap = 2" << std::endl
0438              << "    station = int(station)" << std::endl
0439              << "    ring = int(ring)" << std::endl
0440              << "    chamber = int(chamber)" << std::endl
0441              << "    return endcap, station, ring, chamber" << std::endl
0442              << std::endl;
0443     }
0444 
0445     for (std::vector<CSCChamberFitter>::const_iterator fitter = m_fitters.begin(); fitter != m_fitters.end();
0446          ++fitter) {
0447       if (m_mode == CSCPairResidualsConstraint::kModeRadius) {
0448         fitter->radiusCorrection(m_alignableNavigator, m_alignmentParameterStore, m_combineME11);
0449 
0450       } else {
0451         std::vector<CSCAlignmentCorrections*> corrections;
0452         fitter->fit(corrections);
0453 
0454         // corrections only exist if the fit was successful
0455         for (std::vector<CSCAlignmentCorrections*>::iterator correction = corrections.begin();
0456              correction != corrections.end();
0457              ++correction) {
0458           (*correction)->applyAlignment(m_alignableNavigator, m_alignmentParameterStore, m_mode, m_combineME11);
0459           if (m_makeHistograms)
0460             (*correction)->plot();
0461           if (writeReport)
0462             (*correction)->report(report);
0463         }
0464       }
0465     }
0466   }
0467 }
0468 
0469 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
0470 DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory, CSCOverlapsAlignmentAlgorithm, "CSCOverlapsAlignmentAlgorithm");