00001
00002
00003
00004
00005
00013
00014
00015
00016
00017
00018 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmBase.h"
00019 #include "FWCore/Framework/interface/ESHandle.h"
00020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00022 #include "Alignment/CommonAlignment/interface/Alignable.h"
00023 #include "Alignment/CommonAlignment/interface/AlignableDetUnit.h"
00024 #include "Alignment/MuonAlignment/interface/AlignableDTSuperLayer.h"
00025 #include "Alignment/MuonAlignment/interface/AlignableDTChamber.h"
00026 #include "Alignment/MuonAlignment/interface/AlignableDTStation.h"
00027 #include "Alignment/MuonAlignment/interface/AlignableDTWheel.h"
00028 #include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
00029 #include "Alignment/MuonAlignment/interface/AlignableCSCRing.h"
00030 #include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"
00031 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
00032 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
00033 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
00034 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00035 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00036 #include "DataFormats/MuonDetId/interface/DTSuperLayerId.h"
00037 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
00038 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00039 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00040 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00041 #include "DataFormats/TrackReco/interface/Track.h"
00042 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00043 #include "FWCore/ServiceRegistry/interface/Service.h"
00044 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00045 #include "TFile.h"
00046
00047 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
00048 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
00049 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals5DOFFitter.h"
00050 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFrphiFitter.h"
00051 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h"
00052
00053 #include <map>
00054 #include <sstream>
00055 #include <fstream>
00056
00057 class MuonAlignmentFromReference : public AlignmentAlgorithmBase {
00058 public:
00059 MuonAlignmentFromReference(const edm::ParameterSet& iConfig);
00060 ~MuonAlignmentFromReference();
00061
00062 void initialize(const edm::EventSetup& iSetup,
00063 AlignableTracker* alignableTracker, AlignableMuon* alignableMuon,
00064 AlignableExtras* extras,
00065 AlignmentParameterStore* alignmentParameterStore);
00066 void startNewLoop();
00067 void run(const edm::EventSetup& iSetup, const EventInfo &eventInfo);
00068
00069 void terminate();
00070
00071 private:
00072 bool numeric(std::string s);
00073 int number(std::string s);
00074
00075 std::vector<std::string> m_reference;
00076 double m_minTrackPt;
00077 double m_maxTrackPt;
00078 int m_minTrackerHits;
00079 double m_maxTrackerRedChi2;
00080 bool m_allowTIDTEC;
00081 int m_minDT13Hits;
00082 int m_minDT2Hits;
00083 int m_minCSCHits;
00084 std::string m_writeTemporaryFile;
00085 std::vector<std::string> m_readTemporaryFiles;
00086 bool m_doAlignment;
00087 int m_strategy;
00088 std::string m_residualsModel;
00089 int m_minAlignmentHits;
00090 bool m_twoBin;
00091 bool m_combineME11;
00092 bool m_weightAlignment;
00093 std::string m_reportFileName;
00094 double m_maxResSlopeY;
00095
00096 AlignableNavigator *m_alignableNavigator;
00097 AlignmentParameterStore *m_alignmentParameterStore;
00098 std::vector<Alignable*> m_alignables;
00099 std::map<Alignable*,Alignable*> m_me11map;
00100 std::map<Alignable*,MuonResidualsTwoBin*> m_fitters;
00101 std::vector<unsigned int> m_indexes;
00102 std::map<unsigned int,MuonResidualsTwoBin*> m_fitterOrder;
00103
00104 long m_counter_events;
00105 long m_counter_tracks;
00106 long m_counter_trackpt;
00107 long m_counter_trackerhits;
00108 long m_counter_trackerchi2;
00109 long m_counter_trackertidtec;
00110 long m_counter_station123;
00111 long m_counter_station123valid;
00112 long m_counter_station123dt13hits;
00113 long m_counter_station123dt2hits;
00114 long m_counter_station123aligning;
00115 long m_counter_station4;
00116 long m_counter_station4valid;
00117 long m_counter_station4hits;
00118 long m_counter_station4aligning;
00119 long m_counter_csc;
00120 long m_counter_cscvalid;
00121 long m_counter_cschits;
00122 long m_counter_cscaligning;
00123 long m_counter_resslopey;
00124 };
00125
00126 MuonAlignmentFromReference::MuonAlignmentFromReference(const edm::ParameterSet &iConfig)
00127 : AlignmentAlgorithmBase(iConfig)
00128 , m_reference(iConfig.getParameter<std::vector<std::string> >("reference"))
00129 , m_minTrackPt(iConfig.getParameter<double>("minTrackPt"))
00130 , m_maxTrackPt(iConfig.getParameter<double>("maxTrackPt"))
00131 , m_minTrackerHits(iConfig.getParameter<int>("minTrackerHits"))
00132 , m_maxTrackerRedChi2(iConfig.getParameter<double>("maxTrackerRedChi2"))
00133 , m_allowTIDTEC(iConfig.getParameter<bool>("allowTIDTEC"))
00134 , m_minDT13Hits(iConfig.getParameter<int>("minDT13Hits"))
00135 , m_minDT2Hits(iConfig.getParameter<int>("minDT2Hits"))
00136 , m_minCSCHits(iConfig.getParameter<int>("minCSCHits"))
00137 , m_writeTemporaryFile(iConfig.getParameter<std::string>("writeTemporaryFile"))
00138 , m_readTemporaryFiles(iConfig.getParameter<std::vector<std::string> >("readTemporaryFiles"))
00139 , m_doAlignment(iConfig.getParameter<bool>("doAlignment"))
00140 , m_strategy(iConfig.getParameter<int>("strategy"))
00141 , m_residualsModel(iConfig.getParameter<std::string>("residualsModel"))
00142 , m_minAlignmentHits(iConfig.getParameter<int>("minAlignmentHits"))
00143 , m_twoBin(iConfig.getParameter<bool>("twoBin"))
00144 , m_combineME11(iConfig.getParameter<bool>("combineME11"))
00145 , m_weightAlignment(iConfig.getParameter<bool>("weightAlignment"))
00146 , m_reportFileName(iConfig.getParameter<std::string>("reportFileName"))
00147 , m_maxResSlopeY(iConfig.getParameter<double>("maxResSlopeY"))
00148 {
00149
00150
00151
00152 if (m_doAlignment) {
00153 edm::Service<TFileService> tfileService;
00154 TFile &tfile = tfileService->file();
00155 tfile.ls();
00156 }
00157
00158 m_counter_events = 0;
00159 m_counter_tracks = 0;
00160 m_counter_trackpt = 0;
00161 m_counter_trackerhits = 0;
00162 m_counter_trackerchi2 = 0;
00163 m_counter_trackertidtec = 0;
00164 m_counter_station123 = 0;
00165 m_counter_station123valid = 0;
00166 m_counter_station123dt13hits = 0;
00167 m_counter_station123dt2hits = 0;
00168 m_counter_station123aligning = 0;
00169 m_counter_station4 = 0;
00170 m_counter_station4valid = 0;
00171 m_counter_station4hits = 0;
00172 m_counter_station4aligning = 0;
00173 m_counter_csc = 0;
00174 m_counter_cscvalid = 0;
00175 m_counter_cschits = 0;
00176 m_counter_cscaligning = 0;
00177 m_counter_resslopey = 0;
00178 }
00179
00180 MuonAlignmentFromReference::~MuonAlignmentFromReference() {
00181 delete m_alignableNavigator;
00182 }
00183
00184 bool MuonAlignmentFromReference::numeric(std::string s) {
00185 return (s == std::string("0") || s == std::string("1") || s == std::string("2") || s == std::string("3") || s == std::string("4") ||
00186 s == std::string("5") || s == std::string("6") || s == std::string("7") || s == std::string("8") || s == std::string("9"));
00187 }
00188
00189 int MuonAlignmentFromReference::number(std::string s) {
00190 if (s == std::string("0")) return 0;
00191 else if (s == std::string("1")) return 1;
00192 else if (s == std::string("2")) return 2;
00193 else if (s == std::string("3")) return 3;
00194 else if (s == std::string("4")) return 4;
00195 else if (s == std::string("5")) return 5;
00196 else if (s == std::string("6")) return 6;
00197 else if (s == std::string("7")) return 7;
00198 else if (s == std::string("8")) return 8;
00199 else if (s == std::string("9")) return 9;
00200 else assert(false);
00201 }
00202
00203 void MuonAlignmentFromReference::initialize(const edm::EventSetup& iSetup,
00204 AlignableTracker* alignableTracker, AlignableMuon* alignableMuon,
00205 AlignableExtras* extras,
00206 AlignmentParameterStore* alignmentParameterStore) {
00207 if (alignableMuon == NULL) {
00208 throw cms::Exception("MuonAlignmentFromReference") << "doMuon must be set to True" << std::endl;
00209 }
00210
00211 m_alignableNavigator = new AlignableNavigator(alignableMuon);
00212 m_alignmentParameterStore = alignmentParameterStore;
00213 m_alignables = m_alignmentParameterStore->alignables();
00214
00215 int residualsModel;
00216 if (m_residualsModel == std::string("pureGaussian")) residualsModel = MuonResidualsFitter::kPureGaussian;
00217 else if (m_residualsModel == std::string("powerLawTails")) residualsModel = MuonResidualsFitter::kPowerLawTails;
00218 else if (m_residualsModel == std::string("ROOTVoigt")) residualsModel = MuonResidualsFitter::kROOTVoigt;
00219 else if (m_residualsModel == std::string("GaussPowerTails")) residualsModel = MuonResidualsFitter::kGaussPowerTails;
00220 else throw cms::Exception("MuonAlignmentFromReference") << "unrecognized residualsModel: \"" << m_residualsModel << "\"" << std::endl;
00221
00222 edm::ESHandle<CSCGeometry> cscGeometry;
00223 iSetup.get<MuonGeometryRecord>().get(cscGeometry);
00224
00225
00226 m_me11map.clear();
00227 m_fitters.clear();
00228 m_indexes.clear();
00229 m_fitterOrder.clear();
00230 for (std::vector<Alignable*>::const_iterator ali = m_alignables.begin(); ali != m_alignables.end(); ++ali) {
00231 bool made_fitter = false;
00232
00233 if ((*ali)->alignableObjectId() == align::AlignableDTChamber) {
00234 DTChamberId id((*ali)->geomDetId().rawId());
00235 if (id.station() == 4) {
00236 m_fitters[*ali] = new MuonResidualsTwoBin(m_twoBin, new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, m_weightAlignment), new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, m_weightAlignment));
00237 made_fitter = true;
00238 }
00239 else {
00240 m_fitters[*ali] = new MuonResidualsTwoBin(m_twoBin, new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, m_weightAlignment), new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, m_weightAlignment));
00241 made_fitter = true;
00242 }
00243 }
00244
00245 else if ((*ali)->alignableObjectId() == align::AlignableCSCChamber) {
00246 Alignable *thisali = *ali;
00247 CSCDetId id((*ali)->geomDetId().rawId());
00248
00249 if (m_combineME11 && id.station() == 1 && id.ring() == 4) {
00250 CSCDetId pairid(id.endcap(), 1, 1, id.chamber());
00251
00252 for (std::vector<Alignable*>::const_iterator ali2 = m_alignables.begin(); ali2 != m_alignables.end(); ++ali2) {
00253 if ((*ali2)->alignableObjectId() == align::AlignableCSCChamber && (*ali2)->geomDetId().rawId() == pairid.rawId()) {
00254 thisali = *ali2;
00255 break;
00256 }
00257 }
00258 m_me11map[*ali] = thisali;
00259 }
00260
00261 if (thisali == *ali) {
00262 m_fitters[*ali] = new MuonResidualsTwoBin(m_twoBin, new MuonResiduals6DOFrphiFitter(residualsModel, m_minAlignmentHits, &(*cscGeometry), m_weightAlignment), new MuonResiduals6DOFrphiFitter(residualsModel, m_minAlignmentHits, &(*cscGeometry), m_weightAlignment));
00263 made_fitter = true;
00264 }
00265 }
00266
00267 else {
00268 throw cms::Exception("MuonAlignmentFromReference") << "only DTChambers and CSCChambers can be aligned with this module" << std::endl;
00269 }
00270
00271 if (made_fitter) {
00272 m_fitters[*ali]->setStrategy(m_strategy);
00273
00274 int index = (*ali)->geomDetId().rawId();
00275 m_indexes.push_back(index);
00276 m_fitterOrder[index] = m_fitters[*ali];
00277 }
00278 }
00279
00280
00281 std::sort(m_indexes.begin(), m_indexes.end());
00282
00283
00284 std::vector<Alignable*> all_DT_chambers = alignableMuon->DTChambers();
00285 std::vector<Alignable*> all_CSC_chambers = alignableMuon->CSCChambers();
00286 std::vector<Alignable*> reference;
00287 std::map<Alignable*,bool> already_seen;
00288
00289 for (std::vector<std::string>::const_iterator name = m_reference.begin(); name != m_reference.end(); ++name) {
00290 bool parsing_error = false;
00291
00292 bool barrel = (name->substr(0, 2) == std::string("MB"));
00293 bool endcap = (name->substr(0, 2) == std::string("ME"));
00294 if (!barrel && !endcap) parsing_error = true;
00295
00296 if (!parsing_error && barrel) {
00297 int index = 2;
00298
00299 if (name->substr(index, 1) == std::string(" ")) {
00300 index++;
00301 }
00302
00303 bool plus = true;
00304 if (name->substr(index, 1) == std::string("+")) {
00305 plus = true;
00306 index++;
00307 }
00308 else if (name->substr(index, 1) == std::string("-")) {
00309 plus = false;
00310 index++;
00311 }
00312 else if (numeric(name->substr(index, 1))) {}
00313 else parsing_error = true;
00314
00315 int wheel = 0;
00316 bool wheel_digit = false;
00317 while (!parsing_error && numeric(name->substr(index, 1))) {
00318 wheel *= 10;
00319 wheel += number(name->substr(index, 1));
00320 wheel_digit = true;
00321 index++;
00322 }
00323 if (!plus) wheel *= -1;
00324 if (!wheel_digit) parsing_error = true;
00325
00326 if (name->substr(index, 1) != std::string(" ")) parsing_error = true;
00327 index++;
00328
00329 int station = 0;
00330 bool station_digit = false;
00331 while (!parsing_error && numeric(name->substr(index, 1))) {
00332 station *= 10;
00333 station += number(name->substr(index, 1));
00334 station_digit = true;
00335 index++;
00336 }
00337 if (!station_digit) parsing_error = true;
00338
00339 if (name->substr(index, 1) != std::string(" ")) parsing_error = true;
00340 index++;
00341
00342 int sector = 0;
00343 bool sector_digit = false;
00344 while (!parsing_error && numeric(name->substr(index, 1))) {
00345 sector *= 10;
00346 sector += number(name->substr(index, 1));
00347 sector_digit = true;
00348 index++;
00349 }
00350 if (!sector_digit) parsing_error = true;
00351
00352 if (!parsing_error) {
00353 bool no_such_chamber = false;
00354
00355 if (wheel < -2 || wheel > 2) no_such_chamber = true;
00356 if (station < 1 || station > 4) no_such_chamber = true;
00357 if (station == 4 && (sector < 1 || sector > 14)) no_such_chamber = true;
00358 if (station < 4 && (sector < 1 || sector > 12)) no_such_chamber = true;
00359
00360 if (no_such_chamber) {
00361 throw cms::Exception("MuonAlignmentFromReference") << "reference chamber doesn't exist: " << (*name) << std::endl;
00362 }
00363
00364 DTChamberId id(wheel, station, sector);
00365 for (std::vector<Alignable*>::const_iterator ali = all_DT_chambers.begin(); ali != all_DT_chambers.end(); ++ali) {
00366 if ((*ali)->geomDetId().rawId() == id.rawId()) {
00367 std::map<Alignable*,bool>::const_iterator trial = already_seen.find(*ali);
00368 if (trial == already_seen.end()) {
00369 reference.push_back(*ali);
00370 already_seen[*ali] = true;
00371 }
00372 }
00373 }
00374 }
00375 }
00376 if (!parsing_error && endcap) {
00377 int index = 2;
00378
00379 if (name->substr(index, 1) == std::string(" ")) {
00380 index++;
00381 }
00382
00383 bool plus = true;
00384 if (name->substr(index, 1) == std::string("+")) {
00385 plus = true;
00386 index++;
00387 }
00388 else if (name->substr(index, 1) == std::string("-")) {
00389 plus = false;
00390 index++;
00391 }
00392 else if (numeric(name->substr(index, 1))) {}
00393 else parsing_error = true;
00394
00395 int station = 0;
00396 bool station_digit = false;
00397 while (!parsing_error && numeric(name->substr(index, 1))) {
00398 station *= 10;
00399 station += number(name->substr(index, 1));
00400 station_digit = true;
00401 index++;
00402 }
00403 if (!plus) station *= -1;
00404 if (!station_digit) parsing_error = true;
00405
00406 if (name->substr(index, 1) != std::string("/")) parsing_error = true;
00407 index++;
00408
00409 int ring = 0;
00410 bool ring_digit = false;
00411 while (!parsing_error && numeric(name->substr(index, 1))) {
00412 ring *= 10;
00413 ring += number(name->substr(index, 1));
00414 ring_digit = true;
00415 index++;
00416 }
00417 if (!ring_digit) parsing_error = true;
00418
00419 if (name->substr(index, 1) != std::string(" ")) parsing_error = true;
00420 index++;
00421
00422 int chamber = 0;
00423 bool chamber_digit = false;
00424 while (!parsing_error && numeric(name->substr(index, 1))) {
00425 chamber *= 10;
00426 chamber += number(name->substr(index, 1));
00427 chamber_digit = true;
00428 index++;
00429 }
00430 if (!chamber_digit) parsing_error = true;
00431
00432 if (!parsing_error) {
00433 bool no_such_chamber = false;
00434
00435 int endcap = (station > 0 ? 1 : 2);
00436 station = abs(station);
00437 if (station < 1 || station > 4) no_such_chamber = true;
00438 if (station == 1 && (ring < 1 || ring > 4)) no_such_chamber = true;
00439 if (station > 1 && (ring < 1 || ring > 2)) no_such_chamber = true;
00440 if (station == 1 && (chamber < 1 || chamber > 36)) no_such_chamber = true;
00441 if (station > 1 && ring == 1 && (chamber < 1 || chamber > 18)) no_such_chamber = true;
00442 if (station > 1 && ring == 2 && (chamber < 1 || chamber > 36)) no_such_chamber = true;
00443
00444 if (no_such_chamber) {
00445 throw cms::Exception("MuonAlignmentFromReference") << "reference chamber doesn't exist: " << (*name) << std::endl;
00446 }
00447
00448 CSCDetId id(endcap, station, ring, chamber);
00449 for (std::vector<Alignable*>::const_iterator ali = all_CSC_chambers.begin(); ali != all_CSC_chambers.end(); ++ali) {
00450 if ((*ali)->geomDetId().rawId() == id.rawId()) {
00451 std::map<Alignable*,bool>::const_iterator trial = already_seen.find(*ali);
00452 if (trial == already_seen.end()) {
00453 reference.push_back(*ali);
00454 already_seen[*ali] = true;
00455 }
00456 }
00457 }
00458 }
00459 }
00460
00461 if (parsing_error) {
00462 throw cms::Exception("MuonAlignmentFromReference") << "reference chamber name is malformed: " << (*name) << std::endl;
00463 }
00464 }
00465
00466 alignmentParameterStore->setAlignmentPositionError(all_DT_chambers, 1000., 0.);
00467 alignmentParameterStore->setAlignmentPositionError(all_CSC_chambers, 1000., 0.);
00468 alignmentParameterStore->setAlignmentPositionError(reference, 0., 0.);
00469 }
00470
00471 void MuonAlignmentFromReference::startNewLoop() {}
00472
00473 void MuonAlignmentFromReference::run(const edm::EventSetup& iSetup, const EventInfo &eventInfo) {
00474 m_counter_events++;
00475
00476 edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
00477 iSetup.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
00478
00479 const ConstTrajTrackPairCollection &trajtracks = eventInfo.trajTrackPairs_;
00480 for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end(); ++trajtrack) {
00481 m_counter_tracks++;
00482
00483 const Trajectory* traj = (*trajtrack).first;
00484 const reco::Track* track = (*trajtrack).second;
00485
00486 if (m_minTrackPt < track->pt() && track->pt() < m_maxTrackPt) {
00487 m_counter_trackpt++;
00488
00489 char charge = (track->charge() > 0 ? 1 : -1);
00490
00491
00492 MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, traj, m_alignableNavigator, 1000.);
00493
00494 if (muonResidualsFromTrack.trackerNumHits() >= m_minTrackerHits) {
00495 m_counter_trackerhits++;
00496 if (muonResidualsFromTrack.trackerRedChi2() < m_maxTrackerRedChi2) {
00497 m_counter_trackerchi2++;
00498 if (m_allowTIDTEC || !muonResidualsFromTrack.contains_TIDTEC()) {
00499 m_counter_trackertidtec++;
00500
00501 std::vector<DetId> chamberIds = muonResidualsFromTrack.chamberIds();
00502
00503 for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId) {
00504 if (chamberId->det() == DetId::Muon && chamberId->subdetId() == MuonSubdetId::DT && DTChamberId(chamberId->rawId()).station() != 4) {
00505 MuonChamberResidual *dt13 = muonResidualsFromTrack.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
00506 MuonChamberResidual *dt2 = muonResidualsFromTrack.chamberResidual(*chamberId, MuonChamberResidual::kDT2);
00507
00508 m_counter_station123++;
00509 if (dt13 != NULL && dt2 != NULL) {
00510 m_counter_station123valid++;
00511 if (dt13->numHits() >= m_minDT13Hits) {
00512 m_counter_station123dt13hits++;
00513 if (dt2->numHits() >= m_minDT2Hits) {
00514 m_counter_station123dt2hits++;
00515 std::map<Alignable*,MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(dt13->chamberAlignable());
00516 if (fitter != m_fitters.end()) {
00517 m_counter_station123aligning++;
00518 if (fabs(dt2->resslope()) < m_maxResSlopeY) {
00519 m_counter_resslopey++;
00520
00521 double *residdata = new double[MuonResiduals6DOFFitter::kNData];
00522 residdata[MuonResiduals6DOFFitter::kResidX] = dt13->residual();
00523 residdata[MuonResiduals6DOFFitter::kResidY] = dt2->residual();
00524 residdata[MuonResiduals6DOFFitter::kResSlopeX] = dt13->resslope();
00525 residdata[MuonResiduals6DOFFitter::kResSlopeY] = dt2->resslope();
00526 residdata[MuonResiduals6DOFFitter::kPositionX] = dt13->trackx();
00527 residdata[MuonResiduals6DOFFitter::kPositionY] = dt13->tracky();
00528 residdata[MuonResiduals6DOFFitter::kAngleX] = dt13->trackdxdz();
00529 residdata[MuonResiduals6DOFFitter::kAngleY] = dt13->trackdydz();
00530 residdata[MuonResiduals6DOFFitter::kRedChi2] = (dt13->chi2() + dt2->chi2()) / double(dt13->ndof() + dt2->ndof());
00531 fitter->second->fill(charge, residdata);
00532
00533 }
00534 }
00535 }}}
00536 }
00537
00538 else if (chamberId->det() == DetId::Muon && chamberId->subdetId() == MuonSubdetId::DT && DTChamberId(chamberId->rawId()).station() == 4) {
00539 MuonChamberResidual *dt13 = muonResidualsFromTrack.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
00540
00541 m_counter_station4++;
00542 if (dt13 != NULL) {
00543 m_counter_station4valid++;
00544 if (dt13->numHits() >= m_minDT13Hits) {
00545 m_counter_station4hits++;
00546
00547 std::map<Alignable*,MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(dt13->chamberAlignable());
00548 if (fitter != m_fitters.end()) {
00549 m_counter_station4aligning++;
00550
00551 double *residdata = new double[MuonResiduals5DOFFitter::kNData];
00552 residdata[MuonResiduals5DOFFitter::kResid] = dt13->residual();
00553 residdata[MuonResiduals5DOFFitter::kResSlope] = dt13->resslope();
00554 residdata[MuonResiduals5DOFFitter::kPositionX] = dt13->trackx();
00555 residdata[MuonResiduals5DOFFitter::kPositionY] = dt13->tracky();
00556 residdata[MuonResiduals5DOFFitter::kAngleX] = dt13->trackdxdz();
00557 residdata[MuonResiduals5DOFFitter::kAngleY] = dt13->trackdydz();
00558 residdata[MuonResiduals5DOFFitter::kRedChi2] = dt13->chi2() / double(dt13->ndof());
00559 fitter->second->fill(charge, residdata);
00560
00561 }
00562 }}
00563 }
00564
00565 else if (chamberId->det() == DetId::Muon && chamberId->subdetId() == MuonSubdetId::CSC) {
00566 MuonChamberResidual *csc = muonResidualsFromTrack.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
00567
00568 m_counter_csc++;
00569 if (csc != NULL) {
00570 m_counter_cscvalid++;
00571 if (csc->numHits() >= m_minCSCHits) {
00572 m_counter_cschits++;
00573
00574 Alignable *ali = csc->chamberAlignable();
00575 CSCDetId id(ali->geomDetId().rawId());
00576 if (m_combineME11 && id.station() == 1 && id.ring() == 4) {
00577 ali = m_me11map[ali];
00578 }
00579
00580 std::map<Alignable*,MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(ali);
00581 if (fitter != m_fitters.end()) {
00582 m_counter_cscaligning++;
00583
00584 double *residdata = new double[MuonResiduals6DOFrphiFitter::kNData];
00585 residdata[MuonResiduals6DOFrphiFitter::kResid] = csc->residual();
00586 residdata[MuonResiduals6DOFrphiFitter::kResSlope] = csc->resslope();
00587 residdata[MuonResiduals6DOFrphiFitter::kPositionX] = csc->trackx();
00588 residdata[MuonResiduals6DOFrphiFitter::kPositionY] = csc->tracky();
00589 residdata[MuonResiduals6DOFrphiFitter::kAngleX] = csc->trackdxdz();
00590 residdata[MuonResiduals6DOFrphiFitter::kAngleY] = csc->trackdydz();
00591 residdata[MuonResiduals6DOFrphiFitter::kRedChi2] = csc->chi2() / double(csc->ndof());
00592 fitter->second->fill(charge, residdata);
00593
00594 }
00595 }}
00596 }
00597 else { assert(false); }
00598
00599 }
00600 }}}
00601 }
00602 }
00603 }
00604
00605 void MuonAlignmentFromReference::terminate() {
00606
00607 std::cout << "Counters:" << std::endl
00608 << "COUNT{ events: " << m_counter_events << " }" << std::endl
00609 << "COUNT{ tracks: " << m_counter_tracks << " }" << std::endl
00610 << "COUNT{ trackpt: " << m_counter_trackpt << " }" << std::endl
00611 << "COUNT{ trackerhits: " << m_counter_trackerhits << " }" << std::endl
00612 << "COUNT{ trackerchi2: " << m_counter_trackerchi2 << " }" << std::endl
00613 << "COUNT{ trackertidtec: " << m_counter_trackertidtec << " }" << std::endl
00614 << "COUNT{ station123: " << m_counter_station123 << " }" << std::endl
00615 << "COUNT{ station123valid: " << m_counter_station123valid << " }" << std::endl
00616 << "COUNT{ station123dt13hits: " << m_counter_station123dt13hits << " }" << std::endl
00617 << "COUNT{ station123dt2hits: " << m_counter_station123dt2hits << " }" << std::endl
00618 << "COUNT{ station123aligning: " << m_counter_station123aligning << " }" << std::endl
00619 << "COUNT{ resslopey: " << m_counter_resslopey << " }" << std::endl
00620 << "COUNT{ station4: " << m_counter_station4 << " }" << std::endl
00621 << "COUNT{ station4valid: " << m_counter_station4valid << " }" << std::endl
00622 << "COUNT{ station4hits: " << m_counter_station4hits << " }" << std::endl
00623 << "COUNT{ station4aligning: " << m_counter_station4aligning << " }" << std::endl
00624 << "COUNT{ csc: " << m_counter_csc << " }" << std::endl
00625 << "COUNT{ cscvalid: " << m_counter_cscvalid << " }" << std::endl
00626 << "COUNT{ cschits: " << m_counter_cschits << " }" << std::endl
00627 << "COUNT{ cscaligning: " << m_counter_cscaligning << " }" << std::endl
00628 << "That's all!" << std::endl;
00629
00630
00631 if (m_readTemporaryFiles.size() != 0) {
00632 for (std::vector<std::string>::const_iterator fileName = m_readTemporaryFiles.begin(); fileName != m_readTemporaryFiles.end(); ++fileName) {
00633 FILE *file;
00634 int size;
00635 file = fopen(fileName->c_str(), "r");
00636 if (file == NULL) {
00637 throw cms::Exception("MuonAlignmentFromReference") << "file \"" << *fileName << " can't be opened (doesn't exist?)" << std::endl;
00638 }
00639
00640 fread(&size, sizeof(int), 1, file);
00641 if (int(m_indexes.size()) != size) throw cms::Exception("MuonAlignmentFromReference") << "file \"" << *fileName << "\" has " << size << " fitters, but this job has " << m_indexes.size() << " fitters (probably corresponds to the wrong alignment job)" << std::endl;
00642
00643 int i = 0;
00644 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i) {
00645 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
00646 unsigned int index_toread;
00647 fread(&index_toread, sizeof(unsigned int), 1, file);
00648 if (*index != index_toread) throw cms::Exception("MuonAlignmentFromReference") << "file \"" << *fileName << "\" has index " << index_toread << " at position " << i << ", but this job is expecting " << *index << " (probably corresponds to the wrong alignment job)" << std::endl;
00649 fitter->read(file, i);
00650 }
00651
00652 fclose(file);
00653 }
00654 }
00655
00656
00657
00658 if (m_doAlignment) {
00659 edm::Service<TFileService> tfileService;
00660 TFileDirectory rootDirectory(tfileService->mkdir("MuonAlignmentFromReference"));
00661
00662 std::ofstream report;
00663 bool writeReport = (m_reportFileName != std::string(""));
00664 if (writeReport) {
00665 report.open(m_reportFileName.c_str());
00666 report << "nan = None; NAN = None" << std::endl;
00667 report << "reports = []" << std::endl;
00668 report << "class ValErr:" << std::endl
00669 << " def __init__(self, value, error, antisym):" << std::endl
00670 << " self.value, self.error, self.antisym = value, error, antisym" << std::endl
00671 << "" << std::endl
00672 << " def __repr__(self):" << std::endl
00673 << " if self.antisym == 0.:" << std::endl
00674 << " return \"%g +- %g\" % (self.value, self.error)" << std::endl
00675 << " else:" << std::endl
00676 << " return \"%g +- %g ~ %g\" % (self.value, self.error, self.antisym)" << std::endl
00677 << "" << std::endl
00678 << "class Report:" << std::endl
00679 << " def __init__(self, chamberId, postal_address, name):" << std::endl
00680 << " self.chamberId, self.postal_address, self.name = chamberId, postal_address, name" << std::endl
00681 << " self.status = \"NOFIT\"" << std::endl
00682 << " self.fittype = None" << std::endl
00683 << "" << std::endl
00684 << " def add_parameters(self, deltax, deltay, deltaz, deltaphix, deltaphiy, deltaphiz, loglikelihood, numsegments, sumofweights, redchi2):" << std::endl
00685 << " self.status = \"PASS\"" << std::endl
00686 << " self.deltax, self.deltay, self.deltaz, self.deltaphix, self.deltaphiy, self.deltaphiz = deltax, deltay, deltaz, deltaphix, deltaphiy, deltaphiz" << std::endl
00687 << " self.loglikelihood, self.numsegments, self.sumofweights, self.redchi2 = loglikelihood, numsegments, sumofweights, redchi2" << std::endl
00688 << "" << std::endl
00689 << " def add_stats(self, median_x, median_y, median_dxdz, median_dydz, mean30_x, mean30_y, mean20_dxdz, mean50_dydz, mean15_x, mean15_y, mean10_dxdz, mean25_dydz, wmean30_x, wmean30_y, wmean20_dxdz, wmean50_dydz, wmean15_x, wmean15_y, wmean10_dxdz, wmean25_dydz, stdev30_x, stdev30_y, stdev20_dxdz, stdev50_dydz, stdev15_x, stdev15_y, stdev10_dxdz, stdev25_dydz):" << std::endl
00690 << " self.median_x, self.median_y, self.median_dxdz, self.median_dydz, self.mean30_x, self.mean30_y, self.mean20_dxdz, self.mean50_dydz, self.mean15_x, self.mean15_y, self.mean10_dxdz, self.mean25_dydz, self.wmean30_x, self.wmean30_y, self.wmean20_dxdz, self.wmean50_dydz, self.wmean15_x, self.wmean15_y, self.wmean10_dxdz, self.wmean25_dydz, self.stdev30_x, self.stdev30_y, self.stdev20_dxdz, self.stdev50_dydz, self.stdev15_x, self.stdev15_y, self.stdev10_dxdz, self.stdev25_dydz = median_x, median_y, median_dxdz, median_dydz, mean30_x, mean30_y, mean20_dxdz, mean50_dydz, mean15_x, mean15_y, mean10_dxdz, mean25_dydz, wmean30_x, wmean30_y, wmean20_dxdz, wmean50_dydz, wmean15_x, wmean15_y, wmean10_dxdz, wmean25_dydz, stdev30_x, stdev30_y, stdev20_dxdz, stdev50_dydz, stdev15_x, stdev15_y, stdev10_dxdz, stdev25_dydz" << std::endl
00691 << "" << std::endl
00692 << " def __repr__(self):" << std::endl
00693 << " return \"<Report %s %s %s>\" % (self.postal_address[0], \" \".join(map(str, self.postal_address[1:])), self.status)" << std::endl
00694 << std::endl;
00695 }
00696
00697 for (std::vector<Alignable*>::const_iterator ali = m_alignables.begin(); ali != m_alignables.end(); ++ali) {
00698 std::vector<bool> selector = (*ali)->alignmentParameters()->selector();
00699 bool align_x = selector[0];
00700 bool align_y = selector[1];
00701 bool align_z = selector[2];
00702 bool align_phix = selector[3];
00703 bool align_phiy = selector[4];
00704 bool align_phiz = selector[5];
00705 int numParams = ((align_x ? 1 : 0) + (align_y ? 1 : 0) + (align_z ? 1 : 0) + (align_phix ? 1 : 0) + (align_phiy ? 1 : 0) + (align_phiz ? 1 : 0));
00706
00707
00708 std::vector<int> paramIndex;
00709 int paramIndex_counter = -1;
00710 if (align_x) paramIndex_counter++;
00711 paramIndex.push_back(paramIndex_counter);
00712 if (align_y) paramIndex_counter++;
00713 paramIndex.push_back(paramIndex_counter);
00714 if (align_z) paramIndex_counter++;
00715 paramIndex.push_back(paramIndex_counter);
00716 if (align_phix) paramIndex_counter++;
00717 paramIndex.push_back(paramIndex_counter);
00718 if (align_phiy) paramIndex_counter++;
00719 paramIndex.push_back(paramIndex_counter);
00720 if (align_phiz) paramIndex_counter++;
00721 paramIndex.push_back(paramIndex_counter);
00722
00723 DetId id = (*ali)->geomDetId();
00724
00725 Alignable *thisali = *ali;
00726 if (m_combineME11 && id.subdetId() == MuonSubdetId::CSC) {
00727 CSCDetId cscid(id.rawId());
00728 if (cscid.station() == 1 && cscid.ring() == 4) {
00729 thisali = m_me11map[*ali];
00730 }
00731 }
00732
00733 std::map<Alignable*,MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(thisali);
00734
00735 std::stringstream name;
00736 if (id.subdetId() == MuonSubdetId::DT) {
00737 DTChamberId chamberId(id.rawId());
00738 name << "MBwh";
00739 if (chamberId.wheel() == -2) name << "A";
00740 else if (chamberId.wheel() == -1) name << "B";
00741 else if (chamberId.wheel() == 0) name << "C";
00742 else if (chamberId.wheel() == +1) name << "D";
00743 else if (chamberId.wheel() == +2) name << "E";
00744 std::string sectoro("0");
00745 if (chamberId.sector() > 9) sectoro = std::string("");
00746 name << "st" << chamberId.station() << "sec" << sectoro << chamberId.sector();
00747
00748 if (writeReport) {
00749 report << "reports.append(Report(" << id.rawId() << ", (\"DT\", " << chamberId.wheel() << ", " << chamberId.station() << ", " << chamberId.sector() << "), \"" << name.str() << "\"))" << std::endl;
00750 }
00751 }
00752 else if (id.subdetId() == MuonSubdetId::CSC) {
00753 CSCDetId chamberId(id.rawId());
00754 std::string chambero("0");
00755 if (chamberId.chamber() > 9) chambero = std::string("");
00756 name << "ME" << (chamberId.endcap() == 1 ? "p" : "m") << abs(chamberId.station()) << chamberId.ring() << "_" << chambero << chamberId.chamber();
00757
00758 if (writeReport) {
00759 report << "reports.append(Report(" << id.rawId() << ", (\"CSC\", " << chamberId.endcap() << ", " << chamberId.station() << ", " << chamberId.ring() << ", " << chamberId.chamber() << "), \"" << name.str() << "\"))" << std::endl;
00760 }
00761 }
00762
00763 if (fitter != m_fitters.end()) {
00764
00765 std::cout << "=============================================================================================" << std::endl;
00766 std::cout << "Fitting " << name.str() << std::endl;
00767
00768 if (writeReport) {
00769 report << "reports[-1].posNum = " << fitter->second->numResidualsPos() << std::endl;
00770 report << "reports[-1].negNum = " << fitter->second->numResidualsNeg() << std::endl;
00771 }
00772
00773 if (fitter->second->type() == MuonResidualsFitter::k5DOF) {
00774 if (!align_x) fitter->second->fix(MuonResiduals5DOFFitter::kAlignX);
00775 if (!align_z) fitter->second->fix(MuonResiduals5DOFFitter::kAlignZ);
00776 if (!align_phix) fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiX);
00777 if (!align_phiy) fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiY);
00778 if (!align_phiz) fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiZ);
00779 }
00780 else if (fitter->second->type() == MuonResidualsFitter::k6DOF) {
00781 if (!align_x) fitter->second->fix(MuonResiduals6DOFFitter::kAlignX);
00782 if (!align_y) fitter->second->fix(MuonResiduals6DOFFitter::kAlignY);
00783 if (!align_z) fitter->second->fix(MuonResiduals6DOFFitter::kAlignZ);
00784 if (!align_phix) fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiX);
00785 if (!align_phiy) fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiY);
00786 if (!align_phiz) fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiZ);
00787 }
00788 else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi) {
00789 if (!align_x) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignX);
00790 if (!align_y) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignY);
00791 if (!align_z) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignZ);
00792 if (!align_phix) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiX);
00793 if (!align_phiy) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiY);
00794 if (!align_phiz) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
00795 }
00796 else { assert(false); }
00797
00798 AlgebraicVector params(numParams);
00799 AlgebraicSymMatrix cov(numParams);
00800
00801 if (fitter->second->numsegments() >= m_minAlignmentHits) {
00802 bool successful_fit = fitter->second->fit(thisali);
00803
00804 double loglikelihood = fitter->second->loglikelihood();
00805 double numsegments = fitter->second->numsegments();
00806 double sumofweights = fitter->second->sumofweights();
00807 double redchi2 = fitter->second->plot(name.str(), &rootDirectory, thisali);
00808
00809 if (fitter->second->type() == MuonResidualsFitter::k5DOF) {
00810 double deltax_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignX);
00811 double deltax_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignX);
00812 double deltax_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignX);
00813
00814 double deltaz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignZ);
00815 double deltaz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignZ);
00816 double deltaz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignZ);
00817
00818 double deltaphix_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiX);
00819 double deltaphix_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiX);
00820 double deltaphix_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiX);
00821
00822 double deltaphiy_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiY);
00823 double deltaphiy_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiY);
00824 double deltaphiy_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiY);
00825
00826 double deltaphiz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiZ);
00827 double deltaphiz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiZ);
00828 double deltaphiz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiZ);
00829
00830 double sigmaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidSigma);
00831 double sigmaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidSigma);
00832 double sigmaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidSigma);
00833
00834 double sigmaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeSigma);
00835 double sigmaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeSigma);
00836 double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeSigma);
00837
00838 double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error, gammaresslope_antisym;
00839 gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error = gammaresslope_antisym = 0.;
00840 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian && fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
00841 gammaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidGamma);
00842 gammaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidGamma);
00843 gammaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidGamma);
00844
00845 gammaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeGamma);
00846 gammaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeGamma);
00847 gammaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeGamma);
00848 }
00849
00850 if (writeReport) {
00851 report << "reports[-1].fittype = \"5DOF\"" << std::endl;
00852 report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", " << deltax_antisym << "), \\" << std::endl
00853 << " None, \\" << std::endl
00854 << " ValErr(" << deltaz_value << ", " << deltaz_error << ", " << deltaz_antisym << "), \\" << std::endl
00855 << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", " << deltaphix_antisym << "), \\" << std::endl
00856 << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", " << deltaphiy_antisym << "), \\" << std::endl
00857 << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", " << deltaphiz_antisym << "), \\" << std::endl
00858 << " " << loglikelihood << ", " << numsegments << ", " << sumofweights << ", " << redchi2 << ")" << std::endl;
00859 report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", " << sigmaresid_antisym << ")" << std::endl;
00860 report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error << ", " << sigmaresslope_antisym << ")" << std::endl;
00861 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian && fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
00862 report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", " << gammaresid_antisym << ")" << std::endl;
00863 report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error << ", " << gammaresslope_antisym << ")" << std::endl;
00864 }
00865
00866 report << "reports[-1].add_stats("
00867 << fitter->second->median(MuonResiduals5DOFFitter::kResid) << ", "
00868 << "None, "
00869 << fitter->second->median(MuonResiduals5DOFFitter::kResSlope) << ", "
00870 << "None, "
00871 << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 30.) << ", "
00872 << "None, "
00873 << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 20.) << ", "
00874 << "None, "
00875 << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 15.) << ", "
00876 << "None, "
00877 << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 10.) << ", "
00878 << "None, "
00879 << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 30.) << ", "
00880 << "None, "
00881 << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 20.) << ", "
00882 << "None, "
00883 << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 15.) << ", "
00884 << "None, "
00885 << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 10.) << ", "
00886 << "None, "
00887 << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 30.) << ", "
00888 << "None, "
00889 << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 20.) << ", "
00890 << "None, "
00891 << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 15.) << ", "
00892 << "None, "
00893 << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 10.) << ", "
00894 << "None)" << std::endl;
00895
00896 std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
00897 namesimple_x << name.str() << "_simple_x";
00898 namesimple_dxdz << name.str() << "_simple_dxdz";
00899 nameweighted_x << name.str() << "_weighted_x";
00900 nameweighted_dxdz << name.str() << "_weighted_dxdz";
00901
00902 fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals5DOFFitter::kResid, 10.);
00903 fitter->second->plotsimple(namesimple_dxdz.str(), &rootDirectory, MuonResiduals5DOFFitter::kResSlope, 1000.);
00904
00905 fitter->second->plotweighted(nameweighted_x.str(), &rootDirectory, MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 10.);
00906 fitter->second->plotweighted(nameweighted_dxdz.str(), &rootDirectory, MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 1000.);
00907 }
00908
00909 if (successful_fit) {
00910 if (align_x) params[paramIndex[0]] = deltax_value;
00911 if (align_z) params[paramIndex[2]] = deltaz_value;
00912 if (align_phix) params[paramIndex[3]] = deltaphix_value;
00913 if (align_phiy) params[paramIndex[4]] = deltaphiy_value;
00914 if (align_phiz) params[paramIndex[5]] = deltaphiz_value;
00915 }
00916 }
00917
00918 else if (fitter->second->type() == MuonResidualsFitter::k6DOF) {
00919 double deltax_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignX);
00920 double deltax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignX);
00921 double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignX);
00922
00923 double deltay_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignY);
00924 double deltay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignY);
00925 double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignY);
00926
00927 double deltaz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignZ);
00928 double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignZ);
00929 double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignZ);
00930
00931 double deltaphix_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiX);
00932 double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiX);
00933 double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiX);
00934
00935 double deltaphiy_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiY);
00936 double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiY);
00937 double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiY);
00938
00939 double deltaphiz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiZ);
00940 double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiZ);
00941 double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiZ);
00942
00943 double sigmax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXSigma);
00944 double sigmax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXSigma);
00945 double sigmax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXSigma);
00946
00947 double sigmay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYSigma);
00948 double sigmay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYSigma);
00949 double sigmay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYSigma);
00950
00951 double sigmadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXSigma);
00952 double sigmadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXSigma);
00953 double sigmadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXSigma);
00954
00955 double sigmadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYSigma);
00956 double sigmadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYSigma);
00957 double sigmadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYSigma);
00958
00959 double gammax_value, gammax_error, gammax_antisym, gammay_value, gammay_error, gammay_antisym, gammadxdz_value, gammadxdz_error, gammadxdz_antisym, gammadydz_value, gammadydz_error, gammadydz_antisym;
00960 gammax_value = gammax_error = gammax_antisym = gammay_value = gammay_error = gammay_antisym = gammadxdz_value = gammadxdz_error = gammadxdz_antisym = gammadydz_value = gammadydz_error = gammadydz_antisym = 0.;
00961 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian && fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
00962 gammax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXGamma);
00963 gammax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXGamma);
00964 gammax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXGamma);
00965
00966 gammay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYGamma);
00967 gammay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYGamma);
00968 gammay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYGamma);
00969
00970 gammadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXGamma);
00971 gammadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXGamma);
00972 gammadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXGamma);
00973
00974 gammadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYGamma);
00975 gammadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYGamma);
00976 gammadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYGamma);
00977 }
00978
00979 if (writeReport) {
00980 report << "reports[-1].fittype = \"6DOF\"" << std::endl;
00981 report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", " << deltax_antisym << "), \\" << std::endl
00982 << " ValErr(" << deltay_value << ", " << deltay_error << ", " << deltay_antisym << "), \\" << std::endl
00983 << " ValErr(" << deltaz_value << ", " << deltaz_error << ", " << deltaz_antisym << "), \\" << std::endl
00984 << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", " << deltaphix_antisym << "), \\" << std::endl
00985 << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", " << deltaphiy_antisym << "), \\" << std::endl
00986 << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", " << deltaphiz_antisym << "), \\" << std::endl
00987 << " " << loglikelihood << ", " << numsegments << ", " << sumofweights << ", " << redchi2 << ")" << std::endl;
00988 report << "reports[-1].sigmax = ValErr(" << sigmax_value << ", " << sigmax_error << ", " << sigmax_antisym << ")" << std::endl;
00989 report << "reports[-1].sigmay = ValErr(" << sigmay_value << ", " << sigmay_error << ", " << sigmay_antisym << ")" << std::endl;
00990 report << "reports[-1].sigmadxdz = ValErr(" << sigmadxdz_value << ", " << sigmadxdz_error << ", " << sigmadxdz_antisym << ")" << std::endl;
00991 report << "reports[-1].sigmadydz = ValErr(" << sigmadydz_value << ", " << sigmadydz_error << ", " << sigmadydz_antisym << ")" << std::endl;
00992 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian && fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
00993 report << "reports[-1].gammax = ValErr(" << gammax_value << ", " << gammax_error << ", " << gammax_antisym << ")" << std::endl;
00994 report << "reports[-1].gammay = ValErr(" << gammay_value << ", " << gammay_error << ", " << gammay_antisym << ")" << std::endl;
00995 report << "reports[-1].gammadxdz = ValErr(" << gammadxdz_value << ", " << gammadxdz_error << ", " << gammadxdz_antisym << ")" << std::endl;
00996 report << "reports[-1].gammadydz = ValErr(" << gammadydz_value << ", " << gammadydz_error << ", " << gammadydz_antisym << ")" << std::endl;
00997 }
00998
00999 report << "reports[-1].add_stats("
01000 << fitter->second->median(MuonResiduals6DOFFitter::kResidX) << ", "
01001 << fitter->second->median(MuonResiduals6DOFFitter::kResidY) << ", "
01002 << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeX) << ", "
01003 << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeY) << ", "
01004 << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
01005 << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
01006 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
01007 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
01008 << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
01009 << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
01010 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
01011 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ", "
01012 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 30.) << ", "
01013 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 30.) << ", "
01014 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 20.) << ", "
01015 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 50.) << ", "
01016 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 15.) << ", "
01017 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 15.) << ", "
01018 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 10.) << ", "
01019 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 25.) << ", "
01020 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
01021 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
01022 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
01023 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
01024 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
01025 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
01026 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
01027 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ")" << std::endl;
01028
01029 std::stringstream namesimple_x, namesimple_y, namesimple_dxdz, namesimple_dydz, nameweighted_x, nameweighted_y, nameweighted_dxdz, nameweighted_dydz;
01030 namesimple_x << name.str() << "_simple_x";
01031 namesimple_y << name.str() << "_simple_y";
01032 namesimple_dxdz << name.str() << "_simple_dxdz";
01033 namesimple_dydz << name.str() << "_simple_dydz";
01034 nameweighted_x << name.str() << "_weighted_x";
01035 nameweighted_y << name.str() << "_weighted_y";
01036 nameweighted_dxdz << name.str() << "_weighted_dxdz";
01037 nameweighted_dydz << name.str() << "_weighted_dydz";
01038
01039 fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidX, 10.);
01040 fitter->second->plotsimple(namesimple_y.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidY, 10.);
01041 fitter->second->plotsimple(namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeX, 1000.);
01042 fitter->second->plotsimple(namesimple_dydz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeY, 1000.);
01043
01044 fitter->second->plotweighted(nameweighted_x.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 10.);
01045 fitter->second->plotweighted(nameweighted_y.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 10.);
01046 fitter->second->plotweighted(nameweighted_dxdz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 1000.);
01047 fitter->second->plotweighted(nameweighted_dydz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 1000.);
01048 }
01049
01050 if (successful_fit) {
01051 if (align_x) params[paramIndex[0]] = deltax_value;
01052 if (align_y) params[paramIndex[1]] = deltay_value;
01053 if (align_z) params[paramIndex[2]] = deltaz_value;
01054 if (align_phix) params[paramIndex[3]] = deltaphix_value;
01055 if (align_phiy) params[paramIndex[4]] = deltaphiy_value;
01056 if (align_phiz) params[paramIndex[5]] = deltaphiz_value;
01057 }
01058 }
01059
01060 else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi) {
01061 double deltax_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignX);
01062 double deltax_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignX);
01063 double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignX);
01064
01065 double deltay_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignY);
01066 double deltay_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignY);
01067 double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignY);
01068
01069 double deltaz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignZ);
01070 double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignZ);
01071 double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignZ);
01072
01073 double deltaphix_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiX);
01074 double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiX);
01075 double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiX);
01076
01077 double deltaphiy_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiY);
01078 double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiY);
01079 double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiY);
01080
01081 double deltaphiz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
01082 double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
01083 double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
01084
01085 double sigmaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidSigma);
01086 double sigmaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidSigma);
01087 double sigmaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidSigma);
01088
01089 double sigmaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
01090 double sigmaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
01091 double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
01092
01093 double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error, gammaresslope_antisym;
01094 gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error = gammaresslope_antisym = 0.;
01095 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian && fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
01096 gammaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidGamma);
01097 gammaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidGamma);
01098 gammaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidGamma);
01099
01100 gammaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
01101 gammaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
01102 gammaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
01103 }
01104
01105 if (writeReport) {
01106 report << "reports[-1].fittype = \"6DOFrphi\"" << std::endl;
01107 report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", " << deltax_antisym << "), \\" << std::endl
01108 << " ValErr(" << deltay_value << ", " << deltay_error << ", " << deltay_antisym << "), \\" << std::endl
01109 << " ValErr(" << deltaz_value << ", " << deltaz_error << ", " << deltaz_antisym << "), \\" << std::endl
01110 << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", " << deltaphix_antisym << "), \\" << std::endl
01111 << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", " << deltaphiy_antisym << "), \\" << std::endl
01112 << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", " << deltaphiz_antisym << "), \\" << std::endl
01113 << " " << loglikelihood << ", " << numsegments << ", " << sumofweights << ", " << redchi2 << ")" << std::endl;
01114 report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", " << sigmaresid_antisym << ")" << std::endl;
01115 report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error << ", " << sigmaresslope_antisym << ")" << std::endl;
01116 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian && fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails) {
01117 report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", " << gammaresid_antisym << ")" << std::endl;
01118 report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error << ", " << gammaresslope_antisym << ")" << std::endl;
01119 }
01120
01121 report << "reports[-1].add_stats("
01122 << fitter->second->median(MuonResiduals6DOFrphiFitter::kResid) << ", "
01123 << "None, "
01124 << fitter->second->median(MuonResiduals6DOFrphiFitter::kResSlope) << ", "
01125 << "None, "
01126 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", "
01127 << "None, "
01128 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", "
01129 << "None, "
01130 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", "
01131 << "None, "
01132 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", "
01133 << "None, "
01134 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 30.) << ", "
01135 << "None, "
01136 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 20.) << ", "
01137 << "None, "
01138 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 15.) << ", "
01139 << "None, "
01140 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 10.) << ", "
01141 << "None, "
01142 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", "
01143 << "None, "
01144 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", "
01145 << "None, "
01146 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", "
01147 << "None, "
01148 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", "
01149 << "None)" << std::endl;
01150
01151 std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
01152 namesimple_x << name.str() << "_simple_x";
01153 namesimple_dxdz << name.str() << "_simple_dxdz";
01154 nameweighted_x << name.str() << "_weighted_x";
01155 nameweighted_dxdz << name.str() << "_weighted_dxdz";
01156
01157 fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResid, 10.);
01158 fitter->second->plotsimple(namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResSlope, 1000.);
01159
01160 fitter->second->plotweighted(nameweighted_x.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 10.);
01161 fitter->second->plotweighted(nameweighted_dxdz.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 1000.);
01162 }
01163
01164 if (successful_fit) {
01165 if (align_x) params[paramIndex[0]] = deltax_value;
01166 if (align_y) params[paramIndex[1]] = deltay_value;
01167 if (align_z) params[paramIndex[2]] = deltaz_value;
01168 if (align_phix) params[paramIndex[3]] = deltaphix_value;
01169 if (align_phiy) params[paramIndex[4]] = deltaphiy_value;
01170 if (align_phiz) params[paramIndex[5]] = deltaphiz_value;
01171 }
01172 }
01173
01174 if (successful_fit) {
01175 std::vector<Alignable*> oneortwo;
01176 oneortwo.push_back(*ali);
01177 if (thisali != *ali) oneortwo.push_back(thisali);
01178 m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 0., 0.);
01179 }
01180 else {
01181 std::cout << "MINUIT fit failed!" << std::endl;
01182 if (writeReport) {
01183 report << "reports[-1].status = \"MINUITFAIL\"" << std::endl;
01184 }
01185
01186 for (int i = 0; i < numParams; i++) {
01187 cov[i][i] = 1000.;
01188 }
01189
01190 std::vector<Alignable*> oneortwo;
01191 oneortwo.push_back(*ali);
01192 if (thisali != *ali) oneortwo.push_back(thisali);
01193 m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
01194 }
01195 }
01196 else {
01197 std::cout << "Too few hits!" << std::endl;
01198 if (writeReport) {
01199 report << "reports[-1].status = \"TOOFEWHITS\"" << std::endl;
01200 }
01201
01202 for (int i = 0; i < numParams; i++) {
01203 cov[i][i] = 1000.;
01204 }
01205
01206 std::vector<Alignable*> oneortwo;
01207 oneortwo.push_back(*ali);
01208 if (thisali != *ali) oneortwo.push_back(thisali);
01209 m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
01210 }
01211
01212 AlignmentParameters *parnew = (*ali)->alignmentParameters()->cloneFromSelected(params, cov);
01213 (*ali)->setAlignmentParameters(parnew);
01214 m_alignmentParameterStore->applyParameters(*ali);
01215 (*ali)->alignmentParameters()->setValid(true);
01216
01217 }
01218
01219 if (writeReport) report << std::endl;
01220
01221 }
01222
01223 if (writeReport) report.close();
01224 }
01225
01226
01227 if (m_writeTemporaryFile != std::string("")) {
01228 FILE *file;
01229 file = fopen(m_writeTemporaryFile.c_str(), "w");
01230 int size = m_indexes.size();
01231 fwrite(&size, sizeof(int), 1, file);
01232
01233 int i = 0;
01234 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i) {
01235 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
01236 unsigned int index_towrite = *index;
01237 fwrite(&index_towrite, sizeof(unsigned int), 1, file);
01238 fitter->write(file, i);
01239 }
01240
01241 fclose(file);
01242 }
01243 }
01244
01245 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
01246 DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory, MuonAlignmentFromReference, "MuonAlignmentFromReference");
1.7.3