00001
00002
00003
00004
00005
00013
00014
00015
00016
00017
00018 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmBase.h"
00019
00020 #include "FWCore/Framework/interface/Event.h"
00021 #include "FWCore/Framework/interface/ESHandle.h"
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00024 #include "FWCore/Utilities/interface/InputTag.h"
00025 #include "FWCore/ServiceRegistry/interface/Service.h"
00026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00027
00028 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
00029 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00030 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00031 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00032
00033 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00034 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00035 #include "DataFormats/MuonDetId/interface/DTSuperLayerId.h"
00036 #include "DataFormats/TrackReco/interface/Track.h"
00037 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00038
00039 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00040
00041 #include "Alignment/CommonAlignment/interface/Alignable.h"
00042 #include "Alignment/CommonAlignment/interface/AlignableDetUnit.h"
00043 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
00044 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
00045 #include "Alignment/MuonAlignment/interface/AlignableDTSuperLayer.h"
00046 #include "Alignment/MuonAlignment/interface/AlignableDTChamber.h"
00047 #include "Alignment/MuonAlignment/interface/AlignableDTStation.h"
00048 #include "Alignment/MuonAlignment/interface/AlignableDTWheel.h"
00049 #include "Alignment/MuonAlignment/interface/AlignableCSCChamber.h"
00050 #include "Alignment/MuonAlignment/interface/AlignableCSCRing.h"
00051 #include "Alignment/MuonAlignment/interface/AlignableCSCStation.h"
00052 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
00053
00054 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
00055 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
00056 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals5DOFFitter.h"
00057 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFrphiFitter.h"
00058 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h"
00059
00060 #include "TFile.h"
00061 #include "TTree.h"
00062 #include "TStopwatch.h"
00063
00064 #include <map>
00065 #include <sstream>
00066 #include <fstream>
00067
00068
00069 class MuonAlignmentFromReference : public AlignmentAlgorithmBase
00070 {
00071 public:
00072
00073 MuonAlignmentFromReference(const edm::ParameterSet& cfg);
00074 virtual ~MuonAlignmentFromReference();
00075
00076 void initialize(const edm::EventSetup& iSetup,
00077 AlignableTracker* alignableTracker,
00078 AlignableMuon* alignableMuon,
00079 AlignableExtras* extras,
00080 AlignmentParameterStore* alignmentParameterStore);
00081
00082 void startNewLoop() {};
00083
00084 void run(const edm::EventSetup& iSetup, const EventInfo &eventInfo);
00085
00086 void processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft);
00087
00088 void terminate(const edm::EventSetup& iSetup);
00089
00090 private:
00091 bool numeric(std::string s);
00092 int number(std::string s);
00093 std::string chamberPrettyNameFromId(unsigned int idx);
00094
00095 void parseReference(std::vector<Alignable*> &reference,
00096 std::vector<Alignable*> &all_DT_chambers,
00097 std::vector<Alignable*> &all_CSC_chambers);
00098
00099 void fitAndAlign();
00100 void readTmpFiles();
00101 void writeTmpFiles();
00102
00103 void selectResidualsPeaks();
00104 void correctBField();
00105 void eraseNotSelectedResiduals();
00106
00107 void fillNtuple();
00108
00109
00110 edm::InputTag m_muonCollectionTag;
00111 std::vector<std::string> m_reference;
00112 double m_minTrackPt;
00113 double m_maxTrackPt;
00114 double m_minTrackP;
00115 double m_maxTrackP;
00116 double m_maxDxy;
00117 int m_minTrackerHits;
00118 double m_maxTrackerRedChi2;
00119 bool m_allowTIDTEC;
00120 int m_minNCrossedChambers;
00121 int m_minDT13Hits;
00122 int m_minDT2Hits;
00123 int m_minCSCHits;
00124 std::string m_writeTemporaryFile;
00125 std::vector<std::string> m_readTemporaryFiles;
00126 bool m_doAlignment;
00127 int m_strategy;
00128 std::string m_residualsModel;
00129 int m_minAlignmentHits;
00130 bool m_twoBin;
00131 bool m_combineME11;
00132 bool m_weightAlignment;
00133 std::string m_reportFileName;
00134 double m_maxResSlopeY;
00135 bool m_createNtuple;
00136 double m_peakNSigma;
00137 int m_BFieldCorrection;
00138 bool m_doDT;
00139 bool m_doCSC;
00140 std::string m_useResiduals;
00141
00142
00143 AlignableNavigator *m_alignableNavigator;
00144 AlignmentParameterStore *m_alignmentParameterStore;
00145 std::vector<Alignable*> m_alignables;
00146 std::map<Alignable*,Alignable*> m_me11map;
00147 std::map<Alignable*,MuonResidualsTwoBin*> m_fitters;
00148 std::vector<unsigned int> m_indexes;
00149 std::map<unsigned int,MuonResidualsTwoBin*> m_fitterOrder;
00150
00151
00152 long m_counter_events;
00153 long m_counter_tracks;
00154 long m_counter_trackmomentum;
00155 long m_counter_trackdxy;
00156 long m_counter_trackerhits;
00157 long m_counter_trackerchi2;
00158 long m_counter_trackertidtec;
00159 long m_counter_minchambers;
00160 long m_counter_totchambers;
00161 long m_counter_station123;
00162 long m_counter_station123valid;
00163 long m_counter_station123dt13hits;
00164 long m_counter_station123dt2hits;
00165 long m_counter_station123aligning;
00166 long m_counter_station4;
00167 long m_counter_station4valid;
00168 long m_counter_station4hits;
00169 long m_counter_station4aligning;
00170 long m_counter_csc;
00171 long m_counter_cscvalid;
00172 long m_counter_cschits;
00173 long m_counter_cscaligning;
00174 long m_counter_resslopey;
00175
00176
00177 void bookNtuple();
00178 TTree * m_ttree;
00179 MuonResidualsFitter::MuonAlignmentTreeRow m_tree_row;
00180 };
00181
00182
00183 MuonAlignmentFromReference::MuonAlignmentFromReference(const edm::ParameterSet &cfg)
00184 : AlignmentAlgorithmBase(cfg)
00185 , m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag"))
00186 , m_reference(cfg.getParameter<std::vector<std::string> >("reference"))
00187 , m_minTrackPt(cfg.getParameter<double>("minTrackPt"))
00188 , m_maxTrackPt(cfg.getParameter<double>("maxTrackPt"))
00189 , m_minTrackP(cfg.getParameter<double>("minTrackP"))
00190 , m_maxTrackP(cfg.getParameter<double>("maxTrackP"))
00191 , m_maxDxy(cfg.getParameter<double>("maxDxy"))
00192 , m_minTrackerHits(cfg.getParameter<int>("minTrackerHits"))
00193 , m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2"))
00194 , m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC"))
00195 , m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers"))
00196 , m_minDT13Hits(cfg.getParameter<int>("minDT13Hits"))
00197 , m_minDT2Hits(cfg.getParameter<int>("minDT2Hits"))
00198 , m_minCSCHits(cfg.getParameter<int>("minCSCHits"))
00199 , m_writeTemporaryFile(cfg.getParameter<std::string>("writeTemporaryFile"))
00200 , m_readTemporaryFiles(cfg.getParameter<std::vector<std::string> >("readTemporaryFiles"))
00201 , m_doAlignment(cfg.getParameter<bool>("doAlignment"))
00202 , m_strategy(cfg.getParameter<int>("strategy"))
00203 , m_residualsModel(cfg.getParameter<std::string>("residualsModel"))
00204 , m_minAlignmentHits(cfg.getParameter<int>("minAlignmentHits"))
00205 , m_twoBin(cfg.getParameter<bool>("twoBin"))
00206 , m_combineME11(cfg.getParameter<bool>("combineME11"))
00207 , m_weightAlignment(cfg.getParameter<bool>("weightAlignment"))
00208 , m_reportFileName(cfg.getParameter<std::string>("reportFileName"))
00209 , m_maxResSlopeY(cfg.getParameter<double>("maxResSlopeY"))
00210 , m_createNtuple(cfg.getParameter<bool>("createNtuple"))
00211 , m_peakNSigma(cfg.getParameter<double>("peakNSigma"))
00212 , m_BFieldCorrection(cfg.getParameter<int>("bFieldCorrection"))
00213 , m_doDT(cfg.getParameter<bool>("doDT"))
00214 , m_doCSC(cfg.getParameter<bool>("doCSC"))
00215 , m_useResiduals(cfg.getParameter<std::string>("useResiduals"))
00216 {
00217
00218
00219
00220 if (m_doAlignment || m_createNtuple) {
00221 edm::Service<TFileService> fs;
00222 TFile &tfile = fs->file();
00223 tfile.ls();
00224 }
00225
00226 m_ttree = NULL;
00227 if (m_createNtuple) bookNtuple();
00228
00229 m_counter_events = 0;
00230 m_counter_tracks = 0;
00231 m_counter_trackmomentum = 0;
00232 m_counter_trackdxy = 0;
00233 m_counter_trackerhits = 0;
00234 m_counter_trackerchi2 = 0;
00235 m_counter_trackertidtec = 0;
00236 m_counter_minchambers = 0;
00237 m_counter_totchambers = 0;
00238 m_counter_station123 = 0;
00239 m_counter_station123valid = 0;
00240 m_counter_station123dt13hits = 0;
00241 m_counter_station123dt2hits = 0;
00242 m_counter_station123aligning = 0;
00243 m_counter_station4 = 0;
00244 m_counter_station4valid = 0;
00245 m_counter_station4hits = 0;
00246 m_counter_station4aligning = 0;
00247 m_counter_csc = 0;
00248 m_counter_cscvalid = 0;
00249 m_counter_cschits = 0;
00250 m_counter_cscaligning = 0;
00251 m_counter_resslopey = 0;
00252 }
00253
00254
00255 MuonAlignmentFromReference::~MuonAlignmentFromReference()
00256 {
00257 delete m_alignableNavigator;
00258 }
00259
00260
00261 void MuonAlignmentFromReference::bookNtuple()
00262 {
00263 edm::Service<TFileService> fs;
00264 m_ttree = fs->make<TTree>("mual_ttree", "mual_ttree");
00265 m_ttree->Branch("is_plus", &m_tree_row.is_plus, "is_plus/O");
00266 m_ttree->Branch("is_dt", &m_tree_row.is_dt, "is_dt/O");
00267 m_ttree->Branch("station", &m_tree_row.station, "station/b");
00268 m_ttree->Branch("ring_wheel", &m_tree_row.ring_wheel, "ring_wheel/B");
00269 m_ttree->Branch("sector", &m_tree_row.sector, "sector/b");
00270 m_ttree->Branch("res_x", &m_tree_row.res_x, "res_x/F");
00271 m_ttree->Branch("res_y", &m_tree_row.res_y, "res_y/F");
00272 m_ttree->Branch("res_slope_x", &m_tree_row.res_slope_x, "res_slope_x/F");
00273 m_ttree->Branch("res_slope_y", &m_tree_row.res_slope_y, "res_slope_y/F");
00274 m_ttree->Branch("pos_x",&m_tree_row.pos_x, "pos_x/F");
00275 m_ttree->Branch("pos_y",&m_tree_row.pos_y, "pos_y/F");
00276 m_ttree->Branch("angle_x",&m_tree_row.angle_x, "angle_x/F");
00277 m_ttree->Branch("angle_y",&m_tree_row.angle_y,"angle_y/F");
00278 m_ttree->Branch("pz",&m_tree_row.pz,"pz/F");
00279 m_ttree->Branch("pt",&m_tree_row.pt,"pt/F");
00280 m_ttree->Branch("q",&m_tree_row.q,"q/B");
00281 m_ttree->Branch("select", &m_tree_row.select, "select/O");
00282
00283
00284 }
00285
00286
00287 bool MuonAlignmentFromReference::numeric(std::string s)
00288 {
00289 return s.length()==1 && std::isdigit(s[0]);
00290 }
00291
00292
00293 int MuonAlignmentFromReference::number(std::string s)
00294 {
00295 if (!numeric(s)) assert(false);
00296 return atoi(s.c_str());
00297 }
00298
00299
00300 void MuonAlignmentFromReference::initialize(const edm::EventSetup& iSetup,
00301 AlignableTracker* alignableTracker,
00302 AlignableMuon* alignableMuon,
00303 AlignableExtras* extras,
00304 AlignmentParameterStore* alignmentParameterStore)
00305 {
00306 if (alignableMuon == NULL)
00307 throw cms::Exception("MuonAlignmentFromReference") << "doMuon must be set to True" << std::endl;
00308
00309 m_alignableNavigator = new AlignableNavigator(alignableMuon);
00310 m_alignmentParameterStore = alignmentParameterStore;
00311 m_alignables = m_alignmentParameterStore->alignables();
00312
00313 int residualsModel;
00314 if (m_residualsModel == std::string("pureGaussian")) residualsModel = MuonResidualsFitter::kPureGaussian;
00315 else if (m_residualsModel == std::string("pureGaussian2D")) residualsModel = MuonResidualsFitter::kPureGaussian2D;
00316 else if (m_residualsModel == std::string("powerLawTails")) residualsModel = MuonResidualsFitter::kPowerLawTails;
00317 else if (m_residualsModel == std::string("ROOTVoigt")) residualsModel = MuonResidualsFitter::kROOTVoigt;
00318 else if (m_residualsModel == std::string("GaussPowerTails")) residualsModel = MuonResidualsFitter::kGaussPowerTails;
00319 else throw cms::Exception("MuonAlignmentFromReference") << "unrecognized residualsModel: \"" << m_residualsModel << "\"" << std::endl;
00320
00321 int useResiduals;
00322 if (m_useResiduals == std::string("1111")) useResiduals = MuonResidualsFitter::k1111;
00323 else if (m_useResiduals == std::string("1110")) useResiduals = MuonResidualsFitter::k1110;
00324 else if (m_useResiduals == std::string("1100")) useResiduals = MuonResidualsFitter::k1100;
00325 else if (m_useResiduals == std::string("1010")) useResiduals = MuonResidualsFitter::k1010;
00326 else if (m_useResiduals == std::string("0010")) useResiduals = MuonResidualsFitter::k0010;
00327 else throw cms::Exception("MuonAlignmentFromReference") << "unrecognized useResiduals: \"" << m_useResiduals << "\"" << std::endl;
00328
00329 edm::ESHandle<CSCGeometry> cscGeometry;
00330 iSetup.get<MuonGeometryRecord>().get(cscGeometry);
00331
00332
00333 m_me11map.clear();
00334 m_fitters.clear();
00335 m_indexes.clear();
00336 m_fitterOrder.clear();
00337
00338 for (std::vector<Alignable*>::const_iterator ali = m_alignables.begin(); ali != m_alignables.end(); ++ali)
00339 {
00340 bool made_fitter = false;
00341
00342
00343 if ((*ali)->alignableObjectId() == align::AlignableDTChamber)
00344 {
00345 DTChamberId id( (*ali)->geomDetId().rawId() );
00346
00347 if (id.station() == 4)
00348 {
00349 m_fitters[*ali] =
00350 new MuonResidualsTwoBin(m_twoBin, new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment),
00351 new MuonResiduals5DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment));
00352 made_fitter = true;
00353 }
00354 else
00355 {
00356 m_fitters[*ali] =
00357 new MuonResidualsTwoBin(m_twoBin, new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment),
00358 new MuonResiduals6DOFFitter(residualsModel, m_minAlignmentHits, useResiduals, m_weightAlignment));
00359 made_fitter = true;
00360 }
00361 }
00362
00363
00364 else if ((*ali)->alignableObjectId() == align::AlignableCSCChamber)
00365 {
00366 Alignable *thisali = *ali;
00367 CSCDetId id( (*ali)->geomDetId().rawId() );
00368
00369
00370 if (m_combineME11 && id.station() == 1 && id.ring() == 4)
00371 {
00372 CSCDetId pairid(id.endcap(), 1, 1, id.chamber());
00373
00374 for (std::vector<Alignable*>::const_iterator ali2 = m_alignables.begin(); ali2 != m_alignables.end(); ++ali2)
00375 {
00376 if ((*ali2)->alignableObjectId() == align::AlignableCSCChamber && (*ali2)->geomDetId().rawId() == pairid.rawId())
00377 {
00378 thisali = *ali2;
00379 break;
00380 }
00381 }
00382 m_me11map[*ali] = thisali;
00383 }
00384
00385 if (thisali == *ali)
00386 {
00387 m_fitters[*ali] =
00388 new MuonResidualsTwoBin(m_twoBin, new MuonResiduals6DOFrphiFitter(residualsModel, m_minAlignmentHits, useResiduals, &(*cscGeometry), m_weightAlignment),
00389 new MuonResiduals6DOFrphiFitter(residualsModel, m_minAlignmentHits, useResiduals, &(*cscGeometry), m_weightAlignment));
00390 made_fitter = true;
00391 }
00392 }
00393
00394 else {
00395 throw cms::Exception("MuonAlignmentFromReference") << "only DTChambers and CSCChambers can be aligned with this module" << std::endl;
00396 }
00397
00398 if (made_fitter) {
00399 m_fitters[*ali]->setStrategy(m_strategy);
00400
00401 int index = (*ali)->geomDetId().rawId();
00402 m_indexes.push_back(index);
00403 m_fitterOrder[index] = m_fitters[*ali];
00404 }
00405 }
00406
00407
00408 std::sort(m_indexes.begin(), m_indexes.end());
00409
00410
00411 std::vector<Alignable*> all_DT_chambers = alignableMuon->DTChambers();
00412 std::vector<Alignable*> all_CSC_chambers = alignableMuon->CSCChambers();
00413 std::vector<Alignable*> reference;
00414 if (m_reference.size()) parseReference(reference, all_DT_chambers, all_CSC_chambers);
00415
00416 alignmentParameterStore->setAlignmentPositionError(all_DT_chambers, 100000000., 0.);
00417 alignmentParameterStore->setAlignmentPositionError(all_CSC_chambers, 100000000., 0.);
00418 alignmentParameterStore->setAlignmentPositionError(reference, 0., 0.);
00419 }
00420
00421
00422 void MuonAlignmentFromReference::parseReference(
00423 std::vector<Alignable*> &reference,
00424 std::vector<Alignable*> &all_DT_chambers,
00425 std::vector<Alignable*> &all_CSC_chambers)
00426 {
00427 std::map<Alignable*,bool> already_seen;
00428
00429 for (std::vector<std::string>::const_iterator name = m_reference.begin(); name != m_reference.end(); ++name)
00430 {
00431 bool parsing_error = false;
00432
00433 bool barrel = (name->substr(0, 2) == std::string("MB"));
00434 bool endcap = (name->substr(0, 2) == std::string("ME"));
00435 if (!barrel && !endcap) parsing_error = true;
00436
00437 if (!parsing_error && barrel)
00438 {
00439 int index = 2;
00440 if (name->substr(index, 1) == std::string(" ")) index++;
00441
00442 bool plus = true;
00443 if (name->substr(index, 1) == std::string("+"))
00444 {
00445 plus = true;
00446 index++;
00447 }
00448 else if (name->substr(index, 1) == std::string("-"))
00449 {
00450 plus = false;
00451 index++;
00452 }
00453 else if (numeric(name->substr(index, 1))) {}
00454 else parsing_error = true;
00455
00456 int wheel = 0;
00457 bool wheel_digit = false;
00458 while (!parsing_error && numeric(name->substr(index, 1)))
00459 {
00460 wheel *= 10;
00461 wheel += number(name->substr(index, 1));
00462 wheel_digit = true;
00463 index++;
00464 }
00465 if (!plus) wheel *= -1;
00466 if (!wheel_digit) parsing_error = true;
00467
00468 if (name->substr(index, 1) != std::string(" ")) parsing_error = true;
00469 index++;
00470
00471 int station = 0;
00472 bool station_digit = false;
00473 while (!parsing_error && numeric(name->substr(index, 1)))
00474 {
00475 station *= 10;
00476 station += number(name->substr(index, 1));
00477 station_digit = true;
00478 index++;
00479 }
00480 if (!station_digit) parsing_error = true;
00481
00482 if (name->substr(index, 1) != std::string(" ")) parsing_error = true;
00483 index++;
00484
00485 int sector = 0;
00486 bool sector_digit = false;
00487 while (!parsing_error && numeric(name->substr(index, 1)))
00488 {
00489 sector *= 10;
00490 sector += number(name->substr(index, 1));
00491 sector_digit = true;
00492 index++;
00493 }
00494 if (!sector_digit) parsing_error = true;
00495
00496 if (!parsing_error)
00497 {
00498 bool no_such_chamber = false;
00499
00500 if (wheel < -2 || wheel > 2) no_such_chamber = true;
00501 if (station < 1 || station > 4) no_such_chamber = true;
00502 if (station == 4 && (sector < 1 || sector > 14)) no_such_chamber = true;
00503 if (station < 4 && (sector < 1 || sector > 12)) no_such_chamber = true;
00504
00505 if (no_such_chamber)
00506 throw cms::Exception("MuonAlignmentFromReference") << "reference chamber doesn't exist: " << (*name) << std::endl;
00507
00508 DTChamberId id(wheel, station, sector);
00509 for (std::vector<Alignable*>::const_iterator ali = all_DT_chambers.begin(); ali != all_DT_chambers.end(); ++ali)
00510 {
00511 if ((*ali)->geomDetId().rawId() == id.rawId())
00512 {
00513 std::map<Alignable*,bool>::const_iterator trial = already_seen.find(*ali);
00514 if (trial == already_seen.end())
00515 {
00516 reference.push_back(*ali);
00517 already_seen[*ali] = true;
00518 }
00519 }
00520 }
00521 }
00522 }
00523
00524 if (!parsing_error && endcap)
00525 {
00526 int index = 2;
00527 if (name->substr(index, 1) == std::string(" ")) index++;
00528
00529 bool plus = true;
00530 if (name->substr(index, 1) == std::string("+"))
00531 {
00532 plus = true;
00533 index++;
00534 }
00535 else if (name->substr(index, 1) == std::string("-"))
00536 {
00537 plus = false;
00538 index++;
00539 }
00540 else if (numeric(name->substr(index, 1))) {}
00541 else parsing_error = true;
00542
00543 int station = 0;
00544 bool station_digit = false;
00545 while (!parsing_error && numeric(name->substr(index, 1)))
00546 {
00547 station *= 10;
00548 station += number(name->substr(index, 1));
00549 station_digit = true;
00550 index++;
00551 }
00552 if (!plus) station *= -1;
00553 if (!station_digit) parsing_error = true;
00554
00555 if (name->substr(index, 1) != std::string("/")) parsing_error = true;
00556 index++;
00557
00558 int ring = 0;
00559 bool ring_digit = false;
00560 while (!parsing_error && numeric(name->substr(index, 1)))
00561 {
00562 ring *= 10;
00563 ring += number(name->substr(index, 1));
00564 ring_digit = true;
00565 index++;
00566 }
00567 if (!ring_digit) parsing_error = true;
00568
00569 if (name->substr(index, 1) != std::string(" ")) parsing_error = true;
00570 index++;
00571
00572 int chamber = 0;
00573 bool chamber_digit = false;
00574 while (!parsing_error && numeric(name->substr(index, 1)))
00575 {
00576 chamber *= 10;
00577 chamber += number(name->substr(index, 1));
00578 chamber_digit = true;
00579 index++;
00580 }
00581 if (!chamber_digit) parsing_error = true;
00582
00583 if (!parsing_error)
00584 {
00585 bool no_such_chamber = false;
00586
00587 int endcap = (station > 0 ? 1 : 2);
00588 station = abs(station);
00589 if (station < 1 || station > 4) no_such_chamber = true;
00590 if (station == 1 && (ring < 1 || ring > 4)) no_such_chamber = true;
00591 if (station > 1 && (ring < 1 || ring > 2)) no_such_chamber = true;
00592 if (station == 1 && (chamber < 1 || chamber > 36)) no_such_chamber = true;
00593 if (station > 1 && ring == 1 && (chamber < 1 || chamber > 18)) no_such_chamber = true;
00594 if (station > 1 && ring == 2 && (chamber < 1 || chamber > 36)) no_such_chamber = true;
00595
00596 if (no_such_chamber)
00597 throw cms::Exception("MuonAlignmentFromReference") << "reference chamber doesn't exist: " << (*name) << std::endl;
00598
00599 CSCDetId id(endcap, station, ring, chamber);
00600 for (std::vector<Alignable*>::const_iterator ali = all_CSC_chambers.begin(); ali != all_CSC_chambers.end(); ++ali)
00601 {
00602 if ((*ali)->geomDetId().rawId() == id.rawId())
00603 {
00604 std::map<Alignable*,bool>::const_iterator trial = already_seen.find(*ali);
00605 if (trial == already_seen.end())
00606 {
00607 reference.push_back(*ali);
00608 already_seen[*ali] = true;
00609 }
00610 }
00611 }
00612 }
00613 }
00614
00615 if (parsing_error)
00616 throw cms::Exception("MuonAlignmentFromReference") << "reference chamber name is malformed: " << (*name) << std::endl;
00617 }
00618 }
00619
00620
00621 void MuonAlignmentFromReference::run(const edm::EventSetup& iSetup, const EventInfo &eventInfo)
00622 {
00623 m_counter_events++;
00624
00625 edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
00626 iSetup.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
00627
00628 if (m_muonCollectionTag.label().empty())
00629 {
00630 const ConstTrajTrackPairCollection &trajtracks = eventInfo.trajTrackPairs_;
00631 for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end(); ++trajtrack)
00632 {
00633 m_counter_tracks++;
00634
00635 const Trajectory* traj = (*trajtrack).first;
00636 const reco::Track* track = (*trajtrack).second;
00637
00638 if (m_minTrackPt < track->pt() && track->pt() < m_maxTrackPt && m_minTrackP < track->p() && track->p() < m_maxTrackP)
00639 {
00640 m_counter_trackmomentum++;
00641
00642 if ( fabs(track->dxy(eventInfo.beamSpot_.position())) < m_maxDxy )
00643 {
00644 m_counter_trackdxy++;
00645
00646 MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, traj, track, m_alignableNavigator, 1000.);
00647
00648 processMuonResidualsFromTrack(muonResidualsFromTrack);
00649 }
00650 }
00651 }
00652 }
00653 else
00654 {
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682 }
00683 }
00684
00685
00686 void MuonAlignmentFromReference::processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft)
00687 {
00688 if (mrft.trackerNumHits() >= m_minTrackerHits)
00689 {
00690 m_counter_trackerhits++;
00691 if (mrft.normalizedChi2() < m_maxTrackerRedChi2)
00692 {
00693 m_counter_trackerchi2++;
00694 if (m_allowTIDTEC || !mrft.contains_TIDTEC())
00695 {
00696 m_counter_trackertidtec++;
00697
00698 std::vector<DetId> chamberIds = mrft.chamberIds();
00699
00700 if ((int)chamberIds.size() >= m_minNCrossedChambers)
00701 {
00702 m_counter_minchambers++;
00703
00704 char charge = (mrft.getTrack()->charge() > 0 ? 1 : -1);
00705
00706 for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId)
00707 {
00708 if (chamberId->det() != DetId::Muon) continue;
00709 m_counter_totchambers++;
00710
00711
00712 if (m_doDT &&
00713 chamberId->subdetId() == MuonSubdetId::DT &&
00714 DTChamberId(chamberId->rawId()).station() != 4)
00715 {
00716 MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
00717 MuonChamberResidual *dt2 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT2);
00718
00719 m_counter_station123++;
00720 if (dt13 != NULL && dt2 != NULL)
00721 {
00722 m_counter_station123valid++;
00723 if (dt13->numHits() >= m_minDT13Hits)
00724 {
00725 m_counter_station123dt13hits++;
00726 if (dt2->numHits() >= m_minDT2Hits)
00727 {
00728 m_counter_station123dt2hits++;
00729 std::map<Alignable*,MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(dt13->chamberAlignable());
00730 if (fitter != m_fitters.end())
00731 {
00732 m_counter_station123aligning++;
00733 if (fabs(dt2->resslope()) < m_maxResSlopeY)
00734 {
00735 m_counter_resslopey++;
00736 double *residdata = new double[MuonResiduals6DOFFitter::kNData];
00737 residdata[MuonResiduals6DOFFitter::kResidX] = dt13->residual();
00738 residdata[MuonResiduals6DOFFitter::kResidY] = dt2->residual();
00739 residdata[MuonResiduals6DOFFitter::kResSlopeX] = dt13->resslope();
00740 residdata[MuonResiduals6DOFFitter::kResSlopeY] = dt2->resslope();
00741 residdata[MuonResiduals6DOFFitter::kPositionX] = dt13->trackx();
00742 residdata[MuonResiduals6DOFFitter::kPositionY] = dt13->tracky();
00743 residdata[MuonResiduals6DOFFitter::kAngleX] = dt13->trackdxdz();
00744 residdata[MuonResiduals6DOFFitter::kAngleY] = dt13->trackdydz();
00745 residdata[MuonResiduals6DOFFitter::kRedChi2] = (dt13->chi2() + dt2->chi2()) / double(dt13->ndof() + dt2->ndof());
00746 residdata[MuonResiduals6DOFFitter::kPz] = mrft.getTrack()->pz();
00747 residdata[MuonResiduals6DOFFitter::kPt] = mrft.getTrack()->pt();
00748 residdata[MuonResiduals6DOFFitter::kCharge] = mrft.getTrack()->charge();
00749 fitter->second->fill(charge, residdata);
00750
00751 }
00752 }
00753 }
00754 }
00755 }
00756 }
00757
00758
00759 else if (m_doDT &&
00760 chamberId->subdetId() == MuonSubdetId::DT &&
00761 DTChamberId(chamberId->rawId()).station() == 4)
00762 {
00763 MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
00764
00765 m_counter_station4++;
00766 if (dt13 != NULL)
00767 {
00768 m_counter_station4valid++;
00769 if (dt13->numHits() >= m_minDT13Hits)
00770 {
00771 m_counter_station4hits++;
00772
00773 std::map<Alignable*,MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(dt13->chamberAlignable());
00774 if (fitter != m_fitters.end())
00775 {
00776 m_counter_station4aligning++;
00777
00778 double *residdata = new double[MuonResiduals5DOFFitter::kNData];
00779 residdata[MuonResiduals5DOFFitter::kResid] = dt13->residual();
00780 residdata[MuonResiduals5DOFFitter::kResSlope] = dt13->resslope();
00781 residdata[MuonResiduals5DOFFitter::kPositionX] = dt13->trackx();
00782 residdata[MuonResiduals5DOFFitter::kPositionY] = dt13->tracky();
00783 residdata[MuonResiduals5DOFFitter::kAngleX] = dt13->trackdxdz();
00784 residdata[MuonResiduals5DOFFitter::kAngleY] = dt13->trackdydz();
00785 residdata[MuonResiduals5DOFFitter::kRedChi2] = dt13->chi2() / double(dt13->ndof());
00786 residdata[MuonResiduals5DOFFitter::kPz] = mrft.getTrack()->pz();
00787 residdata[MuonResiduals5DOFFitter::kPt] = mrft.getTrack()->pt();
00788 residdata[MuonResiduals5DOFFitter::kCharge] = mrft.getTrack()->charge();
00789 fitter->second->fill(charge, residdata);
00790
00791 }
00792 }
00793 }
00794 }
00795
00796
00797 else if (m_doCSC &&
00798 chamberId->subdetId() == MuonSubdetId::CSC)
00799 {
00800 MuonChamberResidual *csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
00801 m_counter_csc++;
00802 if (csc != NULL)
00803 {
00804 m_counter_cscvalid++;
00805 if (csc->numHits() >= m_minCSCHits)
00806 {
00807 m_counter_cschits++;
00808 Alignable *ali = csc->chamberAlignable();
00809
00810 CSCDetId id(ali->geomDetId().rawId());
00811 if (m_combineME11 && id.station() == 1 && id.ring() == 4) ali = m_me11map[ali];
00812
00813 std::map<Alignable*,MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(ali);
00814 if (fitter != m_fitters.end())
00815 {
00816 m_counter_cscaligning++;
00817 double *residdata = new double[MuonResiduals6DOFrphiFitter::kNData];
00818 residdata[MuonResiduals6DOFrphiFitter::kResid] = csc->residual();
00819 residdata[MuonResiduals6DOFrphiFitter::kResSlope] = csc->resslope();
00820 residdata[MuonResiduals6DOFrphiFitter::kPositionX] = csc->trackx();
00821 residdata[MuonResiduals6DOFrphiFitter::kPositionY] = csc->tracky();
00822 residdata[MuonResiduals6DOFrphiFitter::kAngleX] = csc->trackdxdz();
00823 residdata[MuonResiduals6DOFrphiFitter::kAngleY] = csc->trackdydz();
00824 residdata[MuonResiduals6DOFrphiFitter::kRedChi2] = csc->chi2() / double(csc->ndof());
00825 residdata[MuonResiduals6DOFrphiFitter::kPz] = mrft.getTrack()->pz();
00826 residdata[MuonResiduals6DOFrphiFitter::kPt] = mrft.getTrack()->pt();
00827 residdata[MuonResiduals6DOFrphiFitter::kCharge] = mrft.getTrack()->charge();
00828 fitter->second->fill(charge, residdata);
00829
00830 }
00831 }
00832 }
00833 }
00834
00835 else if (m_doDT && m_doCSC) assert(false);
00836
00837 }
00838 }
00839 }
00840 }
00841 }
00842 }
00843
00844
00845 void MuonAlignmentFromReference::terminate(const edm::EventSetup& iSetup)
00846 {
00847
00848 std::cout << "Counters:" << std::endl
00849 << "COUNT{ events: " << m_counter_events << " }" << std::endl
00850 << "COUNT{ tracks: " << m_counter_tracks << " }" << std::endl
00851 << "COUNT{ trackppt: " << m_counter_trackmomentum << " }" << std::endl
00852 << "COUNT{ trackdxy: " << m_counter_trackdxy << " }" << std::endl
00853 << "COUNT{ trackerhits: " << m_counter_trackerhits << " }" << std::endl
00854 << "COUNT{ trackerchi2: " << m_counter_trackerchi2 << " }" << std::endl
00855 << "COUNT{ trackertidtec: " << m_counter_trackertidtec << " }" << std::endl
00856 << "COUNT{ minnchambers: " << m_counter_minchambers << " }" << std::endl
00857 << "COUNT{ totchambers: " << m_counter_totchambers << " }" << std::endl
00858 << "COUNT{ station123: " << m_counter_station123 << " }" << std::endl
00859 << "COUNT{ station123valid: " << m_counter_station123valid << " }" << std::endl
00860 << "COUNT{ station123dt13hits: " << m_counter_station123dt13hits << " }" << std::endl
00861 << "COUNT{ station123dt2hits: " << m_counter_station123dt2hits << " }" << std::endl
00862 << "COUNT{ station123aligning: " << m_counter_station123aligning << " }" << std::endl
00863 << "COUNT{ resslopey: " << m_counter_resslopey << " }" << std::endl
00864 << "COUNT{ station4: " << m_counter_station4 << " }" << std::endl
00865 << "COUNT{ station4valid: " << m_counter_station4valid << " }" << std::endl
00866 << "COUNT{ station4hits: " << m_counter_station4hits << " }" << std::endl
00867 << "COUNT{ station4aligning: " << m_counter_station4aligning << " }" << std::endl
00868 << "COUNT{ csc: " << m_counter_csc << " }" << std::endl
00869 << "COUNT{ cscvalid: " << m_counter_cscvalid << " }" << std::endl
00870 << "COUNT{ cschits: " << m_counter_cschits << " }" << std::endl
00871 << "COUNT{ cscaligning: " << m_counter_cscaligning << " }" << std::endl
00872 << "That's all!" << std::endl;
00873
00874 TStopwatch stop_watch;
00875
00876
00877 if (m_readTemporaryFiles.size() != 0)
00878 {
00879 stop_watch.Start();
00880 readTmpFiles();
00881 std::cout <<"readTmpFiles took "<< stop_watch.CpuTime() << " sec" << std::endl;
00882 stop_watch.Stop();
00883 }
00884
00885
00886 if (m_peakNSigma > 0. && m_doAlignment)
00887 {
00888 stop_watch.Start();
00889 selectResidualsPeaks();
00890 std::cout <<"selectResidualsPeaks took "<< stop_watch.CpuTime() << " sec" << std::endl;
00891 stop_watch.Stop();
00892 }
00893
00894 if (m_BFieldCorrection > 0 && m_doAlignment)
00895 {
00896 stop_watch.Start();
00897 correctBField();
00898 std::cout <<"correctBField took "<< stop_watch.CpuTime() << " sec" << std::endl;
00899 stop_watch.Stop();
00900 }
00901
00902
00903 if (m_createNtuple)
00904 {
00905 stop_watch.Start();
00906 fillNtuple();
00907 std::cout <<"fillNtuple took "<< stop_watch.CpuTime() << " sec" << std::endl;
00908 stop_watch.Stop();
00909 }
00910
00911 if (m_doAlignment)
00912 {
00913 stop_watch.Start();
00914 eraseNotSelectedResiduals();
00915 std::cout <<"eraseNotSelectedResiduals took "<< stop_watch.CpuTime() << " sec" << std::endl;
00916 stop_watch.Stop();
00917 }
00918
00919
00920 if (m_doAlignment)
00921 {
00922 stop_watch.Start();
00923 fitAndAlign();
00924 std::cout <<"fitAndAlign took "<< stop_watch.CpuTime() << " sec" << std::endl;
00925 stop_watch.Stop();
00926 }
00927
00928
00929 if (m_writeTemporaryFile != std::string("")) writeTmpFiles();
00930 }
00931
00932
00933 void MuonAlignmentFromReference::fitAndAlign()
00934 {
00935 edm::Service<TFileService> tfileService;
00936 TFileDirectory rootDirectory(tfileService->mkdir("MuonAlignmentFromReference"));
00937
00938 std::ofstream report;
00939 bool writeReport = (m_reportFileName != std::string(""));
00940 if (writeReport)
00941 {
00942 report.open(m_reportFileName.c_str());
00943 report << "nan = None; NAN = None" << std::endl;
00944 report << "reports = []" << std::endl;
00945 report << "class ValErr:" << std::endl
00946 << " def __init__(self, value, error, antisym):" << std::endl
00947 << " self.value, self.error, self.antisym = value, error, antisym" << std::endl
00948 << "" << std::endl
00949 << " def __repr__(self):" << std::endl
00950 << " if self.antisym == 0.:" << std::endl
00951 << " return \"%g +- %g\" % (self.value, self.error)" << std::endl
00952 << " else:" << std::endl
00953 << " return \"%g +- %g ~ %g\" % (self.value, self.error, self.antisym)" << std::endl
00954 << "" << std::endl
00955 << "class Report:" << std::endl
00956 << " def __init__(self, chamberId, postal_address, name):" << std::endl
00957 << " self.chamberId, self.postal_address, self.name = chamberId, postal_address, name" << std::endl
00958 << " self.status = \"NOFIT\"" << std::endl
00959 << " self.fittype = None" << std::endl
00960 << "" << std::endl
00961 << " def add_parameters(self, deltax, deltay, deltaz, deltaphix, deltaphiy, deltaphiz, loglikelihood, numsegments, sumofweights, redchi2):" << std::endl
00962 << " self.status = \"PASS\"" << std::endl
00963 << " self.deltax, self.deltay, self.deltaz, self.deltaphix, self.deltaphiy, self.deltaphiz = deltax, deltay, deltaz, deltaphix, deltaphiy, deltaphiz" << std::endl
00964 << " self.loglikelihood, self.numsegments, self.sumofweights, self.redchi2 = loglikelihood, numsegments, sumofweights, redchi2" << std::endl
00965 << "" << std::endl
00966 << " 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
00967 << " 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
00968 << "" << std::endl
00969 << " def __repr__(self):" << std::endl
00970 << " return \"<Report %s %s %s>\" % (self.postal_address[0], \" \".join(map(str, self.postal_address[1:])), self.status)"<< std::endl
00971 << std::endl;
00972 }
00973
00974 for (std::vector<Alignable*>::const_iterator ali = m_alignables.begin(); ali != m_alignables.end(); ++ali)
00975 {
00976 std::vector<bool> selector = (*ali)->alignmentParameters()->selector();
00977 bool align_x = selector[0];
00978 bool align_y = selector[1];
00979 bool align_z = selector[2];
00980 bool align_phix = selector[3];
00981 bool align_phiy = selector[4];
00982 bool align_phiz = selector[5];
00983 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));
00984
00985
00986 std::vector<int> paramIndex;
00987 int paramIndex_counter = -1;
00988 if (align_x) paramIndex_counter++;
00989 paramIndex.push_back(paramIndex_counter);
00990 if (align_y) paramIndex_counter++;
00991 paramIndex.push_back(paramIndex_counter);
00992 if (align_z) paramIndex_counter++;
00993 paramIndex.push_back(paramIndex_counter);
00994 if (align_phix) paramIndex_counter++;
00995 paramIndex.push_back(paramIndex_counter);
00996 if (align_phiy) paramIndex_counter++;
00997 paramIndex.push_back(paramIndex_counter);
00998 if (align_phiz) paramIndex_counter++;
00999 paramIndex.push_back(paramIndex_counter);
01000
01001 DetId id = (*ali)->geomDetId();
01002
01003 Alignable *thisali = *ali;
01004 if (m_combineME11 && id.subdetId() == MuonSubdetId::CSC)
01005 {
01006 CSCDetId cscid(id.rawId());
01007 if (cscid.station() == 1 && cscid.ring() == 4) thisali = m_me11map[*ali];
01008 }
01009
01010 char cname[40];
01011 char wheel_label[][2]={"A","B","C","D","E"};
01012
01013 if (id.subdetId() == MuonSubdetId::DT)
01014 {
01015 DTChamberId chamberId(id.rawId());
01016
01017
01018
01019 sprintf(cname, "MBwh%sst%dsec%02d", wheel_label[chamberId.wheel()+2], chamberId.station(), chamberId.sector());
01020 if (writeReport)
01021 {
01022 report << "reports.append(Report(" << id.rawId() << ", (\"DT\", "
01023 << chamberId.wheel() << ", " << chamberId.station() << ", " << chamberId.sector() << "), \"" << cname << "\"))" << std::endl;
01024 }
01025 }
01026 else if (id.subdetId() == MuonSubdetId::CSC)
01027 {
01028 CSCDetId chamberId(id.rawId());
01029 sprintf(cname, "ME%s%d%d_%02d", (chamberId.endcap() == 1 ? "p" : "m"), chamberId.station(), chamberId.ring(), chamberId.chamber());
01030
01031
01032
01033 if (writeReport)
01034 {
01035 report << "reports.append(Report(" << id.rawId() << ", (\"CSC\", "
01036 << chamberId.endcap() << ", " << chamberId.station() << ", " << chamberId.ring() << ", " << chamberId.chamber()
01037 << "), \"" << cname << "\"))" << std::endl;
01038 }
01039 }
01040
01041
01042
01043 std::map<Alignable*, MuonResidualsTwoBin*>::const_iterator fitter = m_fitters.find(thisali);
01044
01045 if (fitter != m_fitters.end())
01046 {
01047
01048
01049 TStopwatch stop_watch;
01050 stop_watch.Start();
01051
01052
01053 std::cout << "=============================================================================================" << std::endl;
01054 std::cout << "Fitting " << cname << std::endl;
01055
01056 if (writeReport)
01057 {
01058 report << "reports[-1].posNum = " << fitter->second->numResidualsPos() << std::endl;
01059 report << "reports[-1].negNum = " << fitter->second->numResidualsNeg() << std::endl;
01060 }
01061
01062 if (fitter->second->type() == MuonResidualsFitter::k5DOF)
01063 {
01064 if (!align_x) fitter->second->fix(MuonResiduals5DOFFitter::kAlignX);
01065 if (!align_z) fitter->second->fix(MuonResiduals5DOFFitter::kAlignZ);
01066 if (!align_phix) fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiX);
01067 if (!align_phiy) fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiY);
01068 if (!align_phiz) fitter->second->fix(MuonResiduals5DOFFitter::kAlignPhiZ);
01069 }
01070 else if (fitter->second->type() == MuonResidualsFitter::k6DOF)
01071 {
01072 if (!align_x) fitter->second->fix(MuonResiduals6DOFFitter::kAlignX);
01073 if (!align_y) fitter->second->fix(MuonResiduals6DOFFitter::kAlignY);
01074 if (!align_z) fitter->second->fix(MuonResiduals6DOFFitter::kAlignZ);
01075 if (!align_phix) fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiX);
01076 if (!align_phiy) fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiY);
01077 if (!align_phiz) fitter->second->fix(MuonResiduals6DOFFitter::kAlignPhiZ);
01078 }
01079 else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi)
01080 {
01081 if (!align_x) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignX);
01082 if (!align_y) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignY);
01083 if (!align_z) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignZ);
01084 if (!align_phix) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiX);
01085 if (!align_phiy) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiY);
01086 if (!align_phiz) fitter->second->fix(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
01087 }
01088 else assert(false);
01089
01090 AlgebraicVector params(numParams);
01091 AlgebraicSymMatrix cov(numParams);
01092
01093 if (fitter->second->numsegments() >= m_minAlignmentHits)
01094 {
01095 bool successful_fit = fitter->second->fit(thisali);
01096
01097 double loglikelihood = fitter->second->loglikelihood();
01098 double numsegments = fitter->second->numsegments();
01099 double sumofweights = fitter->second->sumofweights();
01100 double redchi2 = fitter->second->plot(cname, &rootDirectory, thisali);
01101
01102 if (fitter->second->type() == MuonResidualsFitter::k5DOF)
01103 {
01104 double deltax_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignX);
01105 double deltax_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignX);
01106 double deltax_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignX);
01107
01108 double deltaz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignZ);
01109 double deltaz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignZ);
01110 double deltaz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignZ);
01111
01112 double deltaphix_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiX);
01113 double deltaphix_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiX);
01114 double deltaphix_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiX);
01115
01116 double deltaphiy_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiY);
01117 double deltaphiy_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiY);
01118 double deltaphiy_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiY);
01119
01120 double deltaphiz_value = fitter->second->value(MuonResiduals5DOFFitter::kAlignPhiZ);
01121 double deltaphiz_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kAlignPhiZ);
01122 double deltaphiz_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kAlignPhiZ);
01123
01124 double sigmaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidSigma);
01125 double sigmaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidSigma);
01126 double sigmaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidSigma);
01127
01128 double sigmaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeSigma);
01129 double sigmaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeSigma);
01130 double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeSigma);
01131
01132 double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error, gammaresslope_antisym;
01133 gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error = gammaresslope_antisym = 0.;
01134
01135 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
01136 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
01137 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails)
01138 {
01139 gammaresid_value = fitter->second->value(MuonResiduals5DOFFitter::kResidGamma);
01140 gammaresid_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResidGamma);
01141 gammaresid_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResidGamma);
01142
01143 gammaresslope_value = fitter->second->value(MuonResiduals5DOFFitter::kResSlopeGamma);
01144 gammaresslope_error = fitter->second->errorerror(MuonResiduals5DOFFitter::kResSlopeGamma);
01145 gammaresslope_antisym = fitter->second->antisym(MuonResiduals5DOFFitter::kResSlopeGamma);
01146 }
01147
01148 if (writeReport)
01149 {
01150 report << "reports[-1].fittype = \"5DOF\"" << std::endl;
01151 report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", " << deltax_antisym << "), \\" << std::endl
01152 << " None, \\" << std::endl
01153 << " ValErr(" << deltaz_value << ", " << deltaz_error << ", " << deltaz_antisym << "), \\" << std::endl
01154 << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", " << deltaphix_antisym << "), \\" << std::endl
01155 << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", " << deltaphiy_antisym << "), \\" << std::endl
01156 << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", " << deltaphiz_antisym << "), \\" << std::endl
01157 << " " << loglikelihood << ", " << numsegments << ", " << sumofweights << ", " << redchi2 << ")" << std::endl;
01158 report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", " << sigmaresid_antisym << ")" << std::endl;
01159 report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error << ", " << sigmaresslope_antisym << ")" << std::endl;
01160 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
01161 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
01162 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails)
01163 {
01164 report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", " << gammaresid_antisym << ")" << std::endl;
01165 report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error << ", " << gammaresslope_antisym << ")" << std::endl;
01166 }
01167
01168 report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals5DOFFitter::kResid) << ", " << "None, "
01169 << fitter->second->median(MuonResiduals5DOFFitter::kResSlope) << ", " << "None, "
01170 << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 30.) << ", " << "None, "
01171 << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 20.) << ", " << "None, "
01172 << fitter->second->mean(MuonResiduals5DOFFitter::kResid, 15.) << ", " << "None, "
01173 << fitter->second->mean(MuonResiduals5DOFFitter::kResSlope, 10.) << ", " << "None, "
01174 << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 30.) << ", " << "None, "
01175 << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 20.) << ", " << "None, "
01176 << fitter->second->wmean(MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 15.) << ", " << "None, "
01177 << fitter->second->wmean(MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 10.) << ", " << "None, "
01178 << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 30.) << ", " << "None, "
01179 << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 20.) << ", " << "None, "
01180 << fitter->second->stdev(MuonResiduals5DOFFitter::kResid, 15.) << ", " << "None, "
01181 << fitter->second->stdev(MuonResiduals5DOFFitter::kResSlope, 10.) << ", " << "None)" << std::endl;
01182
01183 std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
01184 namesimple_x << cname << "_simple_x";
01185 namesimple_dxdz << cname << "_simple_dxdz";
01186 nameweighted_x << cname << "_weighted_x";
01187 nameweighted_dxdz << cname << "_weighted_dxdz";
01188
01189 fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals5DOFFitter::kResid, 10.);
01190 fitter->second->plotsimple(namesimple_dxdz.str(), &rootDirectory, MuonResiduals5DOFFitter::kResSlope, 1000.);
01191
01192 fitter->second->plotweighted(nameweighted_x.str(), &rootDirectory, MuonResiduals5DOFFitter::kResid, MuonResiduals5DOFFitter::kRedChi2, 10.);
01193 fitter->second->plotweighted(nameweighted_dxdz.str(), &rootDirectory, MuonResiduals5DOFFitter::kResSlope, MuonResiduals5DOFFitter::kRedChi2, 1000.);
01194 }
01195
01196 if (successful_fit)
01197 {
01198 if (align_x) params[paramIndex[0]] = deltax_value;
01199 if (align_z) params[paramIndex[2]] = deltaz_value;
01200 if (align_phix) params[paramIndex[3]] = deltaphix_value;
01201 if (align_phiy) params[paramIndex[4]] = deltaphiy_value;
01202 if (align_phiz) params[paramIndex[5]] = deltaphiz_value;
01203 }
01204 }
01205
01206 else if (fitter->second->type() == MuonResidualsFitter::k6DOF)
01207 {
01208 double deltax_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignX);
01209 double deltax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignX);
01210 double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignX);
01211
01212 double deltay_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignY);
01213 double deltay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignY);
01214 double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignY);
01215
01216 double deltaz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignZ);
01217 double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignZ);
01218 double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignZ);
01219
01220 double deltaphix_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiX);
01221 double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiX);
01222 double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiX);
01223
01224 double deltaphiy_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiY);
01225 double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiY);
01226 double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiY);
01227
01228 double deltaphiz_value = fitter->second->value(MuonResiduals6DOFFitter::kAlignPhiZ);
01229 double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kAlignPhiZ);
01230 double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kAlignPhiZ);
01231
01232 double sigmax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXSigma);
01233 double sigmax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXSigma);
01234 double sigmax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXSigma);
01235
01236 double sigmay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYSigma);
01237 double sigmay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYSigma);
01238 double sigmay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYSigma);
01239
01240 double sigmadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXSigma);
01241 double sigmadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXSigma);
01242 double sigmadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXSigma);
01243
01244 double sigmadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYSigma);
01245 double sigmadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYSigma);
01246 double sigmadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYSigma);
01247
01248 double gammax_value, gammax_error, gammax_antisym, gammay_value, gammay_error, gammay_antisym,
01249 gammadxdz_value, gammadxdz_error, gammadxdz_antisym, gammadydz_value, gammadydz_error, gammadydz_antisym;
01250 gammax_value = gammax_error = gammax_antisym = gammay_value = gammay_error = gammay_antisym = gammadxdz_value
01251 = gammadxdz_error = gammadxdz_antisym = gammadydz_value = gammadydz_error = gammadydz_antisym = 0.;
01252 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
01253 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
01254 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails)
01255 {
01256 gammax_value = fitter->second->value(MuonResiduals6DOFFitter::kResidXGamma);
01257 gammax_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidXGamma);
01258 gammax_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidXGamma);
01259
01260 gammay_value = fitter->second->value(MuonResiduals6DOFFitter::kResidYGamma);
01261 gammay_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResidYGamma);
01262 gammay_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResidYGamma);
01263
01264 gammadxdz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeXGamma);
01265 gammadxdz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeXGamma);
01266 gammadxdz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeXGamma);
01267
01268 gammadydz_value = fitter->second->value(MuonResiduals6DOFFitter::kResSlopeYGamma);
01269 gammadydz_error = fitter->second->errorerror(MuonResiduals6DOFFitter::kResSlopeYGamma);
01270 gammadydz_antisym = fitter->second->antisym(MuonResiduals6DOFFitter::kResSlopeYGamma);
01271 }
01272
01273 if (writeReport)
01274 {
01275 report << "reports[-1].fittype = \"6DOF\"" << std::endl;
01276 report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", " << deltax_antisym << "), \\" << std::endl
01277 << " ValErr(" << deltay_value << ", " << deltay_error << ", " << deltay_antisym << "), \\" << std::endl
01278 << " ValErr(" << deltaz_value << ", " << deltaz_error << ", " << deltaz_antisym << "), \\" << std::endl
01279 << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", " << deltaphix_antisym << "), \\" << std::endl
01280 << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", " << deltaphiy_antisym << "), \\" << std::endl
01281 << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", " << deltaphiz_antisym << "), \\" << std::endl
01282 << " " << loglikelihood << ", " << numsegments << ", " << sumofweights << ", " << redchi2 << ")" << std::endl;
01283 report << "reports[-1].sigmax = ValErr(" << sigmax_value << ", " << sigmax_error << ", " << sigmax_antisym<< ")" << std::endl;
01284 report << "reports[-1].sigmay = ValErr(" << sigmay_value << ", " << sigmay_error << ", " << sigmay_antisym<< ")" << std::endl;
01285 report << "reports[-1].sigmadxdz = ValErr(" << sigmadxdz_value << ", " << sigmadxdz_error << ", "<< sigmadxdz_antisym << ")" << std::endl;
01286 report << "reports[-1].sigmadydz = ValErr(" << sigmadydz_value << ", " << sigmadydz_error << ", "<< sigmadydz_antisym << ")" << std::endl;
01287 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
01288 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
01289 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails)
01290 {
01291 report << "reports[-1].gammax = ValErr(" << gammax_value << ", " << gammax_error << ", " << gammax_antisym << ")" << std::endl;
01292 report << "reports[-1].gammay = ValErr(" << gammay_value << ", " << gammay_error << ", " << gammay_antisym << ")" << std::endl;
01293 report << "reports[-1].gammadxdz = ValErr(" << gammadxdz_value << ", " << gammadxdz_error << ", " << gammadxdz_antisym << ")" << std::endl;
01294 report << "reports[-1].gammadydz = ValErr(" << gammadydz_value << ", " << gammadydz_error << ", " << gammadydz_antisym << ")" << std::endl;
01295 }
01296
01297 report << "reports[-1].add_stats("
01298 << fitter->second->median(MuonResiduals6DOFFitter::kResidX) << ", "
01299 << fitter->second->median(MuonResiduals6DOFFitter::kResidY) << ", "
01300 << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeX) << ", "
01301 << fitter->second->median(MuonResiduals6DOFFitter::kResSlopeY) << ", "
01302 << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
01303 << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
01304 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
01305 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
01306 << fitter->second->mean(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
01307 << fitter->second->mean(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
01308 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
01309 << fitter->second->mean(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ", "
01310 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 30.) << ", "
01311 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 30.) << ", "
01312 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX,MuonResiduals6DOFFitter::kRedChi2, 20.) << ", "
01313 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 50.) << ", "
01314 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 15.) << ", "
01315 << fitter->second->wmean(MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 15.) << ", "
01316 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 10.) << ", "
01317 << fitter->second->wmean(MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 25.) << ", "
01318 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 30.) << ", "
01319 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 30.) << ", "
01320 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 20.) << ", "
01321 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 50.) << ", "
01322 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidX, 15.) << ", "
01323 << fitter->second->stdev(MuonResiduals6DOFFitter::kResidY, 15.) << ", "
01324 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeX, 10.) << ", "
01325 << fitter->second->stdev(MuonResiduals6DOFFitter::kResSlopeY, 25.) << ")" << std::endl;
01326
01327 std::stringstream namesimple_x, namesimple_y, namesimple_dxdz, namesimple_dydz, nameweighted_x,
01328 nameweighted_y, nameweighted_dxdz, nameweighted_dydz;
01329 namesimple_x << cname << "_simple_x";
01330 namesimple_y << cname << "_simple_y";
01331 namesimple_dxdz << cname << "_simple_dxdz";
01332 namesimple_dydz << cname << "_simple_dydz";
01333 nameweighted_x << cname << "_weighted_x";
01334 nameweighted_y << cname << "_weighted_y";
01335 nameweighted_dxdz << cname << "_weighted_dxdz";
01336 nameweighted_dydz << cname << "_weighted_dydz";
01337
01338 fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidX, 10.);
01339 fitter->second->plotsimple(namesimple_y.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidY, 10.);
01340 fitter->second->plotsimple(namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeX, 1000.);
01341 fitter->second->plotsimple(namesimple_dydz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeY,1000.);
01342
01343 fitter->second->plotweighted(nameweighted_x.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kRedChi2, 10.);
01344 fitter->second->plotweighted(nameweighted_y.str(), &rootDirectory, MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kRedChi2, 10.);
01345 fitter->second->plotweighted(nameweighted_dxdz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kRedChi2, 1000.);
01346 fitter->second->plotweighted(nameweighted_dydz.str(), &rootDirectory, MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kRedChi2, 1000.);
01347 }
01348
01349 if (successful_fit)
01350 {
01351 if (align_x) params[paramIndex[0]] = deltax_value;
01352 if (align_y) params[paramIndex[1]] = deltay_value;
01353 if (align_z) params[paramIndex[2]] = deltaz_value;
01354 if (align_phix) params[paramIndex[3]] = deltaphix_value;
01355 if (align_phiy) params[paramIndex[4]] = deltaphiy_value;
01356 if (align_phiz) params[paramIndex[5]] = deltaphiz_value;
01357 }
01358 }
01359
01360 else if (fitter->second->type() == MuonResidualsFitter::k6DOFrphi)
01361 {
01362 double deltax_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignX);
01363 double deltax_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignX);
01364 double deltax_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignX);
01365
01366 double deltay_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignY);
01367 double deltay_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignY);
01368 double deltay_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignY);
01369
01370 double deltaz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignZ);
01371 double deltaz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignZ);
01372 double deltaz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignZ);
01373
01374 double deltaphix_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiX);
01375 double deltaphix_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiX);
01376 double deltaphix_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiX);
01377
01378 double deltaphiy_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiY);
01379 double deltaphiy_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiY);
01380 double deltaphiy_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiY);
01381
01382 double deltaphiz_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
01383 double deltaphiz_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
01384 double deltaphiz_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kAlignPhiZ);
01385
01386 double sigmaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidSigma);
01387 double sigmaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidSigma);
01388 double sigmaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidSigma);
01389
01390 double sigmaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
01391 double sigmaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
01392 double sigmaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeSigma);
01393
01394 double gammaresid_value, gammaresid_error, gammaresid_antisym, gammaresslope_value, gammaresslope_error, gammaresslope_antisym;
01395 gammaresid_value = gammaresid_error = gammaresid_antisym = gammaresslope_value = gammaresslope_error = gammaresslope_antisym = 0.;
01396 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
01397 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
01398 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails)
01399 {
01400 gammaresid_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResidGamma);
01401 gammaresid_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResidGamma);
01402 gammaresid_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResidGamma);
01403
01404 gammaresslope_value = fitter->second->value(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
01405 gammaresslope_error = fitter->second->errorerror(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
01406 gammaresslope_antisym = fitter->second->antisym(MuonResiduals6DOFrphiFitter::kResSlopeGamma);
01407 }
01408
01409 if (writeReport)
01410 {
01411 report << "reports[-1].fittype = \"6DOFrphi\"" << std::endl;
01412 report << "reports[-1].add_parameters(ValErr(" << deltax_value << ", " << deltax_error << ", " << deltax_antisym << "), \\" << std::endl
01413 << " ValErr(" << deltay_value << ", " << deltay_error << ", " << deltay_antisym << "), \\" << std::endl
01414 << " ValErr(" << deltaz_value << ", " << deltaz_error << ", " << deltaz_antisym << "), \\" << std::endl
01415 << " ValErr(" << deltaphix_value << ", " << deltaphix_error << ", " << deltaphix_antisym << "), \\" << std::endl
01416 << " ValErr(" << deltaphiy_value << ", " << deltaphiy_error << ", " << deltaphiy_antisym << "), \\" << std::endl
01417 << " ValErr(" << deltaphiz_value << ", " << deltaphiz_error << ", " << deltaphiz_antisym << "), \\" << std::endl
01418 << " " << loglikelihood << ", " << numsegments << ", " << sumofweights << ", " << redchi2 << ")" << std::endl;
01419 report << "reports[-1].sigmaresid = ValErr(" << sigmaresid_value << ", " << sigmaresid_error << ", " << sigmaresid_antisym << ")" << std::endl;
01420 report << "reports[-1].sigmaresslope = ValErr(" << sigmaresslope_value << ", " << sigmaresslope_error << ", " << sigmaresslope_antisym << ")" << std::endl;
01421 if (fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian &&
01422 fitter->second->residualsModel() != MuonResidualsFitter::kPureGaussian2D &&
01423 fitter->second->residualsModel() != MuonResidualsFitter::kGaussPowerTails)
01424 {
01425 report << "reports[-1].gammaresid = ValErr(" << gammaresid_value << ", " << gammaresid_error << ", " << gammaresid_antisym << ")" << std::endl;
01426 report << "reports[-1].gammaresslope = ValErr(" << gammaresslope_value << ", " << gammaresslope_error << ", " << gammaresslope_antisym << ")" << std::endl;
01427 }
01428
01429 report << "reports[-1].add_stats(" << fitter->second->median(MuonResiduals6DOFrphiFitter::kResid) << ", " << "None, "
01430 << fitter->second->median(MuonResiduals6DOFrphiFitter::kResSlope) << ", " << "None, "
01431 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", " << "None, "
01432 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", " << "None, "
01433 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", " << "None, "
01434 << fitter->second->mean(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", " << "None, "
01435 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 30.) << ", " << "None, "
01436 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 20.) << ", " << "None, "
01437 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 15.) << ", " << "None, "
01438 << fitter->second->wmean(MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 10.) << ", " << "None, "
01439 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 30.) << ", " << "None, "
01440 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 20.) << ", " << "None, "
01441 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResid, 15.) << ", " << "None, "
01442 << fitter->second->stdev(MuonResiduals6DOFrphiFitter::kResSlope, 10.) << ", " << "None)" << std::endl;
01443
01444 std::stringstream namesimple_x, namesimple_dxdz, nameweighted_x, nameweighted_dxdz;
01445 namesimple_x << cname << "_simple_x";
01446 namesimple_dxdz << cname << "_simple_dxdz";
01447 nameweighted_x << cname << "_weighted_x";
01448 nameweighted_dxdz << cname << "_weighted_dxdz";
01449
01450 fitter->second->plotsimple(namesimple_x.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResid, 10.);
01451 fitter->second->plotsimple(namesimple_dxdz.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResSlope, 1000.);
01452
01453 fitter->second->plotweighted(nameweighted_x.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResid, MuonResiduals6DOFrphiFitter::kRedChi2, 10.);
01454 fitter->second->plotweighted(nameweighted_dxdz.str(), &rootDirectory, MuonResiduals6DOFrphiFitter::kResSlope, MuonResiduals6DOFrphiFitter::kRedChi2, 1000.);
01455 }
01456
01457 if (successful_fit)
01458 {
01459 if (align_x) params[paramIndex[0]] = deltax_value;
01460 if (align_y) params[paramIndex[1]] = deltay_value;
01461 if (align_z) params[paramIndex[2]] = deltaz_value;
01462 if (align_phix) params[paramIndex[3]] = deltaphix_value;
01463 if (align_phiy) params[paramIndex[4]] = deltaphiy_value;
01464 if (align_phiz) params[paramIndex[5]] = deltaphiz_value;
01465 }
01466 }
01467
01468 if (successful_fit)
01469 {
01470 std::vector<Alignable*> oneortwo;
01471 oneortwo.push_back(*ali);
01472 if (thisali != *ali) oneortwo.push_back(thisali);
01473 m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 0., 0.);
01474 }
01475 else
01476 {
01477 std::cout << "MINUIT fit failed!" << std::endl;
01478 if (writeReport)
01479 {
01480 report << "reports[-1].status = \"MINUITFAIL\"" << std::endl;
01481 }
01482
01483 for (int i = 0; i < numParams; i++) cov[i][i] = 1000.;
01484
01485 std::vector<Alignable*> oneortwo;
01486 oneortwo.push_back(*ali);
01487 if (thisali != *ali) oneortwo.push_back(thisali);
01488 m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
01489 }
01490 }
01491 else
01492 {
01493 std::cout << "Too few hits!" << std::endl;
01494 if (writeReport)
01495 {
01496 report << "reports[-1].status = \"TOOFEWHITS\"" << std::endl;
01497 }
01498
01499 for (int i = 0; i < numParams; i++) cov[i][i] = 1000.;
01500
01501 std::vector<Alignable*> oneortwo;
01502 oneortwo.push_back(*ali);
01503 if (thisali != *ali) oneortwo.push_back(thisali);
01504 m_alignmentParameterStore->setAlignmentPositionError(oneortwo, 1000., 0.);
01505 }
01506
01507 AlignmentParameters *parnew = (*ali)->alignmentParameters()->cloneFromSelected(params, cov);
01508 (*ali)->setAlignmentParameters(parnew);
01509 m_alignmentParameterStore->applyParameters(*ali);
01510 (*ali)->alignmentParameters()->setValid(true);
01511
01512 std::cout << cname<<" fittime= "<< stop_watch.CpuTime() << " sec" << std::endl;
01513 }
01514
01515 if (writeReport) report << std::endl;
01516
01517 }
01518
01519 if (writeReport) report.close();
01520 }
01521
01522
01523 void MuonAlignmentFromReference::readTmpFiles()
01524 {
01525 for (std::vector<std::string>::const_iterator fileName = m_readTemporaryFiles.begin(); fileName != m_readTemporaryFiles.end(); ++fileName)
01526 {
01527 FILE *file;
01528 int size;
01529 file = fopen(fileName->c_str(), "r");
01530 if (file == NULL)
01531 throw cms::Exception("MuonAlignmentFromReference") << "file \"" << *fileName << " can't be opened (doesn't exist?)" << std::endl;
01532
01533 fread(&size, sizeof(int), 1, file);
01534 if (int(m_indexes.size()) != size)
01535 throw cms::Exception("MuonAlignmentFromReference") << "file \"" << *fileName << "\" has " << size
01536 << " fitters, but this job has " << m_indexes.size() << " fitters (probably corresponds to the wrong alignment job)" << std::endl;
01537
01538 int i = 0;
01539 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i)
01540 {
01541 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
01542 unsigned int index_toread;
01543 fread(&index_toread, sizeof(unsigned int), 1, file);
01544 if (*index != index_toread)
01545 throw cms::Exception("MuonAlignmentFromReference") << "file \"" << *fileName << "\" has index " << index_toread
01546 << " at position " << i << ", but this job is expecting " << *index << " (probably corresponds to the wrong alignment job)" << std::endl;
01547 fitter->read(file, i);
01548 }
01549
01550 fclose(file);
01551 }
01552 }
01553
01554
01555 void MuonAlignmentFromReference::writeTmpFiles()
01556 {
01557 FILE *file;
01558 file = fopen(m_writeTemporaryFile.c_str(), "w");
01559 int size = m_indexes.size();
01560 fwrite(&size, sizeof(int), 1, file);
01561
01562 int i = 0;
01563 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index, ++i)
01564 {
01565 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
01566 unsigned int index_towrite = *index;
01567 fwrite(&index_towrite, sizeof(unsigned int), 1, file);
01568 fitter->write(file, i);
01569 }
01570
01571 fclose(file);
01572 }
01573
01574 void MuonAlignmentFromReference::correctBField()
01575 {
01576 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index)
01577 {
01578 std::cout<<"correcting B in "<<chamberPrettyNameFromId(*index)<<std::endl;
01579 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
01580 fitter->correctBField();
01581 }
01582 }
01583
01584
01585 void MuonAlignmentFromReference::eraseNotSelectedResiduals()
01586 {
01587 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index)
01588 {
01589 std::cout<<"erasing in "<<chamberPrettyNameFromId(*index)<<std::endl;
01590 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
01591 fitter->eraseNotSelectedResiduals();
01592 }
01593 }
01594
01595
01596 void MuonAlignmentFromReference::selectResidualsPeaks()
01597 {
01598
01599 assert(m_peakNSigma>0.);
01600
01601 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index)
01602 {
01603 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
01604
01605 int nvar = 2;
01606 int vars_index[10] = {0,1};
01607 if (fitter->type() == MuonResidualsFitter::k5DOF)
01608 {
01609 if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 || fitter->useRes() == MuonResidualsFitter::k1010) {
01610 nvar = 2;
01611 vars_index[0] = MuonResiduals5DOFFitter::kResid;
01612 vars_index[1] = MuonResiduals5DOFFitter::kResSlope;
01613 }
01614 else if (fitter->useRes() == MuonResidualsFitter::k1100) {
01615 nvar = 1;
01616 vars_index[0] = MuonResiduals5DOFFitter::kResid;
01617 }
01618 else if (fitter->useRes() == MuonResidualsFitter::k0010) {
01619 nvar = 1;
01620 vars_index[0] = MuonResiduals5DOFFitter::kResSlope;
01621 }
01622 }
01623 else if (fitter->type() == MuonResidualsFitter::k6DOF)
01624 {
01625 if (fitter->useRes() == MuonResidualsFitter::k1111) {
01626 nvar = 4;
01627 vars_index[0] = MuonResiduals6DOFFitter::kResidX;
01628 vars_index[1] = MuonResiduals6DOFFitter::kResidY;
01629 vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
01630 vars_index[3] = MuonResiduals6DOFFitter::kResSlopeY;
01631 }
01632 else if (fitter->useRes() == MuonResidualsFitter::k1110) {
01633 nvar = 3;
01634 vars_index[0] = MuonResiduals6DOFFitter::kResidX;
01635 vars_index[1] = MuonResiduals6DOFFitter::kResidY;
01636 vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
01637 }
01638 else if (fitter->useRes() == MuonResidualsFitter::k1010) {
01639 nvar = 2;
01640 vars_index[0] = MuonResiduals6DOFFitter::kResidX;
01641 vars_index[2] = MuonResiduals6DOFFitter::kResSlopeX;
01642 }
01643 else if (fitter->useRes() == MuonResidualsFitter::k1100) {
01644 nvar = 2;
01645 vars_index[0] = MuonResiduals6DOFFitter::kResidX;
01646 vars_index[1] = MuonResiduals6DOFFitter::kResidY;
01647 }
01648 else if (fitter->useRes() == MuonResidualsFitter::k0010) {
01649 nvar = 1;
01650 vars_index[0] = MuonResiduals6DOFFitter::kResSlopeX;
01651 }
01652 }
01653 else if (fitter->type() == MuonResidualsFitter::k6DOFrphi)
01654 {
01655 if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 || fitter->useRes() == MuonResidualsFitter::k1010) {
01656 nvar = 2;
01657 vars_index[0] = MuonResiduals6DOFrphiFitter::kResid;
01658 vars_index[1] = MuonResiduals6DOFrphiFitter::kResSlope;
01659 }
01660 else if (fitter->useRes() == MuonResidualsFitter::k1100) {
01661 nvar = 1;
01662 vars_index[0] = MuonResiduals6DOFrphiFitter::kResid;
01663 }
01664 else if (fitter->useRes() == MuonResidualsFitter::k0010) {
01665 nvar = 1;
01666 vars_index[0] = MuonResiduals6DOFrphiFitter::kResSlope;
01667 }
01668 }
01669 else assert(false);
01670
01671 std::cout<<"selecting in "<<chamberPrettyNameFromId(*index)<<std::endl;
01672
01673 fitter->selectPeakResiduals(m_peakNSigma, nvar, vars_index);
01674 }
01675 }
01676
01677
01678 std::string MuonAlignmentFromReference::chamberPrettyNameFromId(unsigned int idx)
01679 {
01680 DetId id(idx);
01681 char cname[40];
01682 if (id.subdetId() == MuonSubdetId::DT)
01683 {
01684 DTChamberId chamberId(id.rawId());
01685 sprintf(cname, "MB%+d/%d/%02d", chamberId.wheel(), chamberId.station(), chamberId.sector());
01686 }
01687 else if (id.subdetId() == MuonSubdetId::CSC)
01688 {
01689 CSCDetId chamberId(id.rawId());
01690 sprintf(cname, "ME%s%d/%d/%02d", (chamberId.endcap() == 1 ? "+" : "-"), chamberId.station(), chamberId.ring(), chamberId.chamber());
01691 }
01692 return std::string(cname);
01693 }
01694
01695
01696 void MuonAlignmentFromReference::fillNtuple()
01697 {
01698
01699
01700 for (std::vector<unsigned int>::const_iterator index = m_indexes.begin(); index != m_indexes.end(); ++index)
01701 {
01702 DetId detid(*index);
01703 if (detid.det() != DetId::Muon || !( detid.subdetId() == MuonSubdetId::DT || detid.subdetId() == MuonSubdetId::CSC) ) assert(false);
01704
01705 if(detid.subdetId() == MuonSubdetId::DT)
01706 {
01707 m_tree_row.is_dt = (Bool_t) true;
01708 DTChamberId id(*index);
01709 m_tree_row.is_plus = (Bool_t) true;
01710 m_tree_row.station = (UChar_t) id.station();
01711 m_tree_row.ring_wheel = (Char_t) id.wheel();
01712 m_tree_row.sector = (UChar_t) id.sector();
01713 }
01714 else
01715 {
01716 m_tree_row.is_dt = (Bool_t) false;
01717 CSCDetId id(*index);
01718 m_tree_row.is_plus = (Bool_t) (id.endcap() == 1);
01719 m_tree_row.station = (UChar_t) id.station();
01720 m_tree_row.ring_wheel = (Char_t) id.ring();
01721 m_tree_row.sector = (UChar_t) id.chamber();
01722 }
01723
01724 MuonResidualsTwoBin *fitter = m_fitterOrder[*index];
01725
01726 std::vector<double*>::const_iterator residual = fitter->residualsPos_begin();
01727 std::vector<bool>::const_iterator residual_ok = fitter->residualsPos_ok_begin();
01728 for (; residual != fitter->residualsPos_end(); ++residual, ++residual_ok)
01729 {
01730 if (fitter->type() == MuonResidualsFitter::k5DOF || fitter->type() == MuonResidualsFitter::k6DOFrphi)
01731 {
01732 m_tree_row.res_x = (Float_t) (*residual)[MuonResiduals5DOFFitter::kResid];
01733 m_tree_row.res_y = (Float_t) 0.;
01734 m_tree_row.res_slope_x = (Float_t) (*residual)[MuonResiduals5DOFFitter::kResSlope];
01735 m_tree_row.res_slope_y = (Float_t) 0.;
01736 m_tree_row.pos_x = (Float_t) (*residual)[MuonResiduals5DOFFitter::kPositionX];
01737 m_tree_row.pos_y = (Float_t) (*residual)[MuonResiduals5DOFFitter::kPositionY];
01738 m_tree_row.angle_x = (Float_t) (*residual)[MuonResiduals5DOFFitter::kAngleX];
01739 m_tree_row.angle_y = (Float_t) (*residual)[MuonResiduals5DOFFitter::kAngleY];
01740 m_tree_row.pz = (Float_t) (*residual)[MuonResiduals5DOFFitter::kPz];
01741 m_tree_row.pt = (Float_t) (*residual)[MuonResiduals5DOFFitter::kPt];
01742 m_tree_row.q = (Char_t) (*residual)[MuonResiduals5DOFFitter::kCharge];
01743 m_tree_row.select = (Bool_t) *residual_ok;
01744 }
01745 else if (fitter->type() == MuonResidualsFitter::k6DOF)
01746 {
01747 m_tree_row.res_x = (Float_t) (*residual)[MuonResiduals6DOFFitter::kResidX];
01748 m_tree_row.res_y = (Float_t) (*residual)[MuonResiduals6DOFFitter::kResidY];
01749 m_tree_row.res_slope_x = (Float_t) (*residual)[MuonResiduals6DOFFitter::kResSlopeX];
01750 m_tree_row.res_slope_y = (Float_t) (*residual)[MuonResiduals6DOFFitter::kResSlopeY];
01751 m_tree_row.pos_x = (Float_t) (*residual)[MuonResiduals6DOFFitter::kPositionX];
01752 m_tree_row.pos_y = (Float_t) (*residual)[MuonResiduals6DOFFitter::kPositionY];
01753 m_tree_row.angle_x = (Float_t) (*residual)[MuonResiduals6DOFFitter::kAngleX];
01754 m_tree_row.angle_y = (Float_t) (*residual)[MuonResiduals6DOFFitter::kAngleY];
01755 m_tree_row.pz = (Float_t) (*residual)[MuonResiduals6DOFFitter::kPz];
01756 m_tree_row.pt = (Float_t) (*residual)[MuonResiduals6DOFFitter::kPt];
01757 m_tree_row.q = (Char_t) (*residual)[MuonResiduals6DOFFitter::kCharge];
01758 m_tree_row.select = (Bool_t) *residual_ok;
01759 }
01760 else assert(false);
01761
01762 m_ttree->Fill();
01763 }
01764 }
01765 }
01766
01767
01768 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
01769 DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory, MuonAlignmentFromReference, "MuonAlignmentFromReference");
1.7.3