File indexing completed on 2024-09-07 04:34:30
0001 #include <fstream>
0002 #include <unordered_map>
0003
0004 #include "TFile.h"
0005 #include "TTree.h"
0006 #include "TList.h"
0007 #include "TRandom.h"
0008 #include "TMath.h"
0009
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Framework/interface/Run.h"
0015 #include "FWCore/Framework/interface/ESTransientHandle.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/Framework/interface/ValidityInterval.h"
0018
0019 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0020 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0021
0022 #include "Alignment/CommonAlignment/interface/Alignable.h"
0023 #include "Alignment/CommonAlignment/interface/AlignableDetUnit.h"
0024 #include "Alignment/CommonAlignment/interface/AlignableExtras.h"
0025 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
0026 #include "Alignment/CommonAlignment/interface/SurveyResidual.h"
0027 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterSelector.h"
0028 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
0029 #include "Alignment/HIPAlignmentAlgorithm/interface/HIPUserVariables.h"
0030 #include "Alignment/HIPAlignmentAlgorithm/interface/HIPUserVariablesIORoot.h"
0031 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
0032 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0033 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0034 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0035 #include "DataFormats/GeometrySurface/interface/LocalError.h"
0036 #include "DataFormats/TrackReco/interface/Track.h"
0037
0038 #include "Alignment/HIPAlignmentAlgorithm/interface/HIPAlignmentAlgorithm.h"
0039
0040
0041 HIPAlignmentAlgorithm::HIPAlignmentAlgorithm(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0042 : AlignmentAlgorithmBase(cfg, iC),
0043 topoToken_(iC.esConsumes<TrackerTopology, IdealGeometryRecord, edm::Transition::EndRun>()),
0044 topoToken2_(iC.esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::EndRun>()),
0045 verbose(cfg.getParameter<bool>("verbosity")),
0046 theMonitorConfig(cfg),
0047 doTrackHitMonitoring(theMonitorConfig.fillTrackMonitoring || theMonitorConfig.fillTrackHitMonitoring),
0048 defaultAlignableSpecs((Alignable*)nullptr),
0049 surveyResiduals_(cfg.getUntrackedParameter<std::vector<std::string>>("surveyResiduals")),
0050 theTrackHitMonitorIORootFile(nullptr),
0051 theTrackMonitorTree(nullptr),
0052 theHitMonitorTree(nullptr),
0053 theAlignablesMonitorIORootFile(nullptr),
0054 theAlignablesMonitorTree(nullptr),
0055 theSurveyIORootFile(nullptr),
0056 theSurveyTree(nullptr) {
0057
0058 outpath = cfg.getParameter<std::string>("outpath");
0059 outfile2 = cfg.getParameter<std::string>("outfile2");
0060 struefile = cfg.getParameter<std::string>("trueFile");
0061 smisalignedfile = cfg.getParameter<std::string>("misalignedFile");
0062 salignedfile = cfg.getParameter<std::string>("alignedFile");
0063 siterationfile = cfg.getParameter<std::string>("iterationFile");
0064 suvarfilecore = cfg.getParameter<std::string>("uvarFile");
0065 suvarfile = suvarfilecore;
0066 sparameterfile = cfg.getParameter<std::string>("parameterFile");
0067 ssurveyfile = cfg.getParameter<std::string>("surveyFile");
0068
0069 outfile2 = outpath + outfile2;
0070 struefile = outpath + struefile;
0071 smisalignedfile = outpath + smisalignedfile;
0072 salignedfile = outpath + salignedfile;
0073 siterationfile = outpath + siterationfile;
0074 suvarfile = outpath + suvarfile;
0075 sparameterfile = outpath + sparameterfile;
0076 ssurveyfile = outpath + ssurveyfile;
0077
0078
0079 theApplyAPE = cfg.getParameter<bool>("applyAPE");
0080 theAPEParameterSet = cfg.getParameter<std::vector<edm::ParameterSet>>("apeParam");
0081
0082 themultiIOV = cfg.getParameter<bool>("multiIOV");
0083 theIOVrangeSet = cfg.getParameter<std::vector<unsigned>>("IOVrange");
0084
0085 defaultAlignableSpecs.minNHits = cfg.getParameter<int>("minimumNumberOfHits");
0086 defaultAlignableSpecs.minRelParError = cfg.getParameter<double>("minRelParameterError");
0087 defaultAlignableSpecs.maxRelParError = cfg.getParameter<double>("maxRelParameterError");
0088 defaultAlignableSpecs.maxHitPull = cfg.getParameter<double>("maxAllowedHitPull");
0089 theApplyCutsPerComponent = cfg.getParameter<bool>("applyCutsPerComponent");
0090 theCutsPerComponent = cfg.getParameter<std::vector<edm::ParameterSet>>("cutsPerComponent");
0091
0092
0093 isCollector = cfg.getParameter<bool>("collectorActive");
0094 theCollectorNJobs = cfg.getParameter<int>("collectorNJobs");
0095 theCollectorPath = cfg.getParameter<std::string>("collectorPath");
0096
0097 if (isCollector)
0098 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::HIPAlignmentAlgorithm"
0099 << "Collector mode";
0100
0101 trackPs = cfg.getParameter<bool>("UsePreSelection");
0102 theDataGroup = cfg.getParameter<int>("DataGroup");
0103 trackWt = cfg.getParameter<bool>("UseReweighting");
0104 Scale = cfg.getParameter<double>("Weight");
0105 uniEta = trackWt && cfg.getParameter<bool>("UniformEta");
0106 uniEtaFormula = cfg.getParameter<std::string>("UniformEtaFormula");
0107 if (uniEtaFormula.empty()) {
0108 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::HIPAlignmentAlgorithm"
0109 << "Uniform eta formula is empty! Resetting to 1.";
0110 uniEtaFormula = "1";
0111 }
0112 theEtaFormula = std::make_unique<TFormula>(uniEtaFormula.c_str());
0113 rewgtPerAli = trackWt && cfg.getParameter<bool>("ReweightPerAlignable");
0114 IsCollision = cfg.getParameter<bool>("isCollision");
0115 SetScanDet = cfg.getParameter<std::vector<double>>("setScanDet");
0116 col_cut = cfg.getParameter<double>("CLAngleCut");
0117 cos_cut = cfg.getParameter<double>("CSAngleCut");
0118
0119 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::HIPAlignmentAlgorithm"
0120 << "Constructed";
0121 }
0122
0123
0124 void HIPAlignmentAlgorithm::initialize(const edm::EventSetup& setup,
0125 AlignableTracker* tracker,
0126 AlignableMuon* muon,
0127 AlignableExtras* extras,
0128 AlignmentParameterStore* store) {
0129 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::initialize"
0130 << "Initializing...";
0131
0132 alignableObjectId_ = std::make_unique<AlignableObjectId>(AlignableObjectId::commonObjectIdProvider(tracker, muon));
0133
0134 for (const auto& level : surveyResiduals_)
0135 theLevels.push_back(alignableObjectId_->stringToId(level));
0136
0137 const edm::ValidityInterval& validity = setup.get<TrackerAlignmentRcd>().validityInterval();
0138 const edm::IOVSyncValue first1 = validity.first();
0139 unsigned int firstrun = first1.eventID().run();
0140 if (themultiIOV) {
0141 if (theIOVrangeSet.size() != 1) {
0142 bool findMatchIOV = false;
0143 for (unsigned int iovl = 0; iovl < theIOVrangeSet.size(); iovl++) {
0144 if (firstrun == theIOVrangeSet.at(iovl)) {
0145 std::string iovapp = std::to_string(firstrun);
0146 iovapp.append(".root");
0147 iovapp.insert(0, "_");
0148 salignedfile.replace(salignedfile.end() - 5, salignedfile.end(), iovapp);
0149 siterationfile.replace(siterationfile.end() - 5, siterationfile.end(), iovapp);
0150
0151 if (isCollector) {
0152 outfile2.replace(outfile2.end() - 5, outfile2.end(), iovapp);
0153 ssurveyfile.replace(ssurveyfile.end() - 5, ssurveyfile.end(), iovapp);
0154 suvarfile.replace(suvarfile.end() - 5, suvarfile.end(), iovapp);
0155 }
0156
0157 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::initialize"
0158 << "Found the IOV file matching IOV first run " << firstrun;
0159 findMatchIOV = true;
0160 break;
0161 }
0162 }
0163 if (!findMatchIOV)
0164 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::initialize"
0165 << "Didn't find the IOV file matching IOV first run " << firstrun
0166 << " from the validity interval";
0167 } else {
0168 std::string iovapp = std::to_string(theIOVrangeSet.at(0));
0169 iovapp.append(".root");
0170 iovapp.insert(0, "_");
0171 salignedfile.replace(salignedfile.end() - 5, salignedfile.end(), iovapp);
0172 siterationfile.replace(siterationfile.end() - 5, siterationfile.end(), iovapp);
0173 }
0174 }
0175
0176
0177 theAlignableDetAccessor = std::make_unique<AlignableNavigator>(extras, tracker, muon);
0178 if (extras != nullptr)
0179 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::initialize"
0180 << "AlignableNavigator initialized with AlignableExtras";
0181
0182
0183 theAlignmentParameterStore = store;
0184
0185
0186 theAlignables = theAlignmentParameterStore->alignables();
0187
0188
0189 {
0190 AlignmentParameterSelector selector(tracker, muon, extras);
0191
0192
0193 theAPEParameters.clear();
0194 if (theApplyAPE) {
0195 for (std::vector<edm::ParameterSet>::const_iterator setiter = theAPEParameterSet.begin();
0196 setiter != theAPEParameterSet.end();
0197 ++setiter) {
0198 align::Alignables alignables;
0199
0200 selector.clear();
0201 edm::ParameterSet selectorPSet = setiter->getParameter<edm::ParameterSet>("Selector");
0202 std::vector<std::string> alignParams = selectorPSet.getParameter<std::vector<std::string>>("alignParams");
0203 if (alignParams.size() == 1 && alignParams[0] == std::string("selected"))
0204 alignables = theAlignables;
0205 else {
0206 selector.addSelections(selectorPSet);
0207 alignables = selector.selectedAlignables();
0208 }
0209
0210 std::vector<double> apeSPar = setiter->getParameter<std::vector<double>>("apeSPar");
0211 std::vector<double> apeRPar = setiter->getParameter<std::vector<double>>("apeRPar");
0212 std::string function = setiter->getParameter<std::string>("function");
0213
0214 if (apeSPar.size() != 3 || apeRPar.size() != 3)
0215 throw cms::Exception("BadConfig") << "apeSPar and apeRPar must have 3 values each" << std::endl;
0216
0217 for (std::vector<double>::const_iterator i = apeRPar.begin(); i != apeRPar.end(); ++i)
0218 apeSPar.push_back(*i);
0219
0220 if (function == std::string("linear"))
0221 apeSPar.push_back(0);
0222 else if (function == std::string("exponential"))
0223 apeSPar.push_back(1);
0224 else if (function == std::string("step"))
0225 apeSPar.push_back(2);
0226 else
0227 throw cms::Exception("BadConfig")
0228 << "APE function must be \"linear\", \"exponential\", or \"step\"." << std::endl;
0229
0230 theAPEParameters.push_back(std::make_pair(alignables, apeSPar));
0231 }
0232 }
0233
0234
0235 theAlignableSpecifics.clear();
0236 if (theApplyCutsPerComponent) {
0237 for (std::vector<edm::ParameterSet>::const_iterator setiter = theCutsPerComponent.begin();
0238 setiter != theCutsPerComponent.end();
0239 ++setiter) {
0240 align::Alignables alignables;
0241
0242 selector.clear();
0243 edm::ParameterSet selectorPSet = setiter->getParameter<edm::ParameterSet>("Selector");
0244 std::vector<std::string> alignParams = selectorPSet.getParameter<std::vector<std::string>>("alignParams");
0245 if (alignParams.size() == 1 && alignParams[0] == std::string("selected"))
0246 alignables = theAlignables;
0247 else {
0248 selector.addSelections(selectorPSet);
0249 alignables = selector.selectedAlignables();
0250 }
0251
0252 double minRelParError = setiter->getParameter<double>("minRelParError");
0253 double maxRelParError = setiter->getParameter<double>("maxRelParError");
0254 int minNHits = setiter->getParameter<int>("minNHits");
0255 double maxHitPull = setiter->getParameter<double>("maxHitPull");
0256 bool applyPixelProbCut = setiter->getParameter<bool>("applyPixelProbCut");
0257 bool usePixelProbXYOrProbQ = setiter->getParameter<bool>("usePixelProbXYOrProbQ");
0258 double minPixelProbXY = setiter->getParameter<double>("minPixelProbXY");
0259 double maxPixelProbXY = setiter->getParameter<double>("maxPixelProbXY");
0260 double minPixelProbQ = setiter->getParameter<double>("minPixelProbQ");
0261 double maxPixelProbQ = setiter->getParameter<double>("maxPixelProbQ");
0262 for (auto& ali : alignables) {
0263 HIPAlignableSpecificParameters alispecs(ali);
0264 alispecs.minRelParError = minRelParError;
0265 alispecs.maxRelParError = maxRelParError;
0266 alispecs.minNHits = minNHits;
0267 alispecs.maxHitPull = maxHitPull;
0268
0269 alispecs.applyPixelProbCut = applyPixelProbCut;
0270 alispecs.usePixelProbXYOrProbQ = usePixelProbXYOrProbQ;
0271 alispecs.minPixelProbXY = minPixelProbXY;
0272 alispecs.maxPixelProbXY = maxPixelProbXY;
0273 alispecs.minPixelProbQ = minPixelProbQ;
0274 alispecs.maxPixelProbQ = maxPixelProbQ;
0275
0276 theAlignableSpecifics.push_back(alispecs);
0277 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::initialize"
0278 << "Alignment specifics acquired for detector " << alispecs.id() << " / "
0279 << alispecs.objId() << ":\n"
0280 << " - minRelParError = " << alispecs.minRelParError << "\n"
0281 << " - maxRelParError = " << alispecs.maxRelParError << "\n"
0282 << " - minNHits = " << alispecs.minNHits << "\n"
0283 << " - maxHitPull = " << alispecs.maxHitPull << "\n"
0284 << " - applyPixelProbCut = " << alispecs.applyPixelProbCut << "\n"
0285 << " - usePixelProbXYOrProbQ = " << alispecs.usePixelProbXYOrProbQ << "\n"
0286 << " - minPixelProbXY = " << alispecs.minPixelProbXY << "\n"
0287 << " - maxPixelProbXY = " << alispecs.maxPixelProbXY << "\n"
0288 << " - minPixelProbQ = " << alispecs.minPixelProbQ << "\n"
0289 << " - maxPixelProbQ = " << alispecs.maxPixelProbQ;
0290 }
0291 }
0292 }
0293 }
0294 }
0295
0296
0297 void HIPAlignmentAlgorithm::startNewLoop(void) {
0298 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0299 << "Begin";
0300
0301
0302 for (const auto& it : theAlignables) {
0303 AlignmentParameters* ap = it->alignmentParameters();
0304 int npar = ap->numSelected();
0305 HIPUserVariables* userpar = new HIPUserVariables(npar);
0306 ap->setUserVariables(userpar);
0307 }
0308
0309
0310 AlignablePositions theAlignablePositionsFromFile =
0311 theIO.readAlignableAbsolutePositions(theAlignables, salignedfile.c_str(), -1, ioerr);
0312 int numAlignablesFromFile = theAlignablePositionsFromFile.size();
0313 if (numAlignablesFromFile == 0) {
0314
0315 if (isCollector)
0316 theIteration = 0;
0317 else
0318 theIteration = 1;
0319 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0320 << "IO alignables file not found for iteration " << theIteration;
0321 } else {
0322 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0323 << "Alignables Read " << numAlignablesFromFile;
0324
0325
0326 theIteration = readIterationFile(siterationfile);
0327
0328 theIO.readAlignableAbsolutePositions(theAlignables, salignedfile.c_str(), theIteration, ioerr);
0329
0330
0331 if (ioerr == 0) {
0332 theIteration++;
0333 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0334 << "Iteration increased by one and is now " << theIteration;
0335 }
0336
0337
0338 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0339 << "Apply positions from file ...";
0340 theAlignmentParameterStore->applyAlignableAbsolutePositions(theAlignables, theAlignablePositionsFromFile, ioerr);
0341 }
0342
0343 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0344 << "Current Iteration number: " << theIteration;
0345
0346
0347 bookRoot();
0348
0349
0350 setAlignmentPositionError();
0351
0352
0353 if (isCollector)
0354 collector();
0355
0356 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0357 << "End";
0358 }
0359
0360
0361 void HIPAlignmentAlgorithm::terminate(const edm::EventSetup& iSetup) {
0362 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Terminating";
0363
0364
0365 if (!theLevels.empty()) {
0366 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint";
0367
0368 unsigned int nAlignable = theAlignables.size();
0369 const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0370 for (unsigned int i = 0; i < nAlignable; ++i) {
0371 const Alignable* ali = theAlignables[i];
0372 AlignmentParameters* ap = ali->alignmentParameters();
0373 HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(ap->userVariables());
0374 int nhit = uservar->nhit;
0375
0376
0377 std::pair<int, int> tl = theAlignmentParameterStore->typeAndLayer(ali, tTopo);
0378 int tmp_Type = tl.first;
0379 int tmp_Layer = tl.second;
0380 GlobalPoint pos = ali->surface().position();
0381 float tmpz = pos.z();
0382 if (nhit < 1500 ||
0383 (tmp_Type == 5 && tmp_Layer == 4 && fabs(tmpz) > 90)) {
0384 for (unsigned int l = 0; l < theLevels.size(); ++l) {
0385 SurveyResidual res(*ali, theLevels[l], true);
0386
0387 if (res.valid()) {
0388 AlgebraicSymMatrix invCov = res.inverseCovariance();
0389
0390
0391 AlgebraicVector sensResid = res.sensorResidual();
0392 m3_Id = ali->id();
0393 m3_ObjId = theLevels[l];
0394 m3_par[0] = sensResid[0];
0395 m3_par[1] = sensResid[1];
0396 m3_par[2] = sensResid[2];
0397 m3_par[3] = sensResid[3];
0398 m3_par[4] = sensResid[4];
0399 m3_par[5] = sensResid[5];
0400
0401 uservar->jtvj += invCov;
0402 uservar->jtve += invCov * sensResid;
0403
0404 if (theSurveyTree != nullptr)
0405 theSurveyTree->Fill();
0406 }
0407 }
0408 }
0409 }
0410 }
0411
0412
0413 HIPUserVariablesIORoot HIPIO;
0414
0415 if (!isCollector)
0416 HIPIO.writeHIPUserVariables(theAlignables, suvarfile.c_str(), theIteration, false, ioerr);
0417
0418
0419 int ialigned = 0;
0420
0421 for (const auto& ali : theAlignables) {
0422 AlignmentParameters* par = ali->alignmentParameters();
0423
0424 if (SetScanDet.at(0) != 0) {
0425 edm::LogWarning("Alignment") << "********Starting Scan*********";
0426 edm::LogWarning("Alignment") << "det ID=" << SetScanDet.at(0) << ", starting position=" << SetScanDet.at(1)
0427 << ", step=" << SetScanDet.at(2) << ", currentDet = " << ali->id();
0428 }
0429
0430 if ((SetScanDet.at(0) != 0) && (SetScanDet.at(0) != 1) && (ali->id() != SetScanDet.at(0)))
0431 continue;
0432
0433 bool test = calcParameters(ali, SetScanDet.at(0), SetScanDet.at(1), SetScanDet.at(2));
0434 if (test) {
0435 if (dynamic_cast<AlignableDetUnit*>(ali) != nullptr) {
0436 std::vector<std::pair<int, SurfaceDeformation*>> pairs;
0437 ali->surfaceDeformationIdPairs(pairs);
0438 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::terminate"
0439 << "The alignable contains " << pairs.size() << " surface deformations";
0440 } else
0441 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::terminate"
0442 << "The alignable cannot contain surface deformations";
0443
0444 theAlignmentParameterStore->applyParameters(ali);
0445
0446 ali->alignmentParameters()->setValid(true);
0447
0448 ialigned++;
0449 } else
0450 par->setValid(false);
0451 }
0452
0453
0454 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::terminate] Aligned units: " << ialigned;
0455
0456
0457 fillAlignablesMonitor(iSetup);
0458
0459 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Writing aligned parameters to file: " << theAlignables.size()
0460 << ", for Iteration " << theIteration;
0461
0462
0463 if (isCollector)
0464 HIPIO.writeHIPUserVariables(theAlignables, suvarfile.c_str(), theIteration, false, ioerr);
0465
0466
0467 theIO.writeAlignableAbsolutePositions(theAlignables, salignedfile.c_str(), theIteration, false, ioerr);
0468
0469
0470
0471
0472
0473
0474 writeIterationFile(siterationfile, theIteration);
0475
0476
0477
0478 if (theSurveyIORootFile != nullptr) {
0479 theSurveyIORootFile->cd();
0480 if (theSurveyTree != nullptr)
0481 theSurveyTree->Write();
0482 delete theSurveyTree;
0483 theSurveyTree = nullptr;
0484 theSurveyIORootFile->Close();
0485 }
0486
0487 if (theAlignablesMonitorIORootFile != nullptr) {
0488 theAlignablesMonitorIORootFile->cd();
0489 if (theAlignablesMonitorTree != nullptr)
0490 theAlignablesMonitorTree->Write();
0491 delete theAlignablesMonitorTree;
0492 theAlignablesMonitorTree = nullptr;
0493 theAlignablesMonitorIORootFile->Close();
0494 }
0495
0496 if (theTrackHitMonitorIORootFile != nullptr) {
0497 theTrackHitMonitorIORootFile->cd();
0498 if (theTrackMonitorTree != nullptr) {
0499 theTrackMonitorTree->Write();
0500 delete theTrackMonitorTree;
0501 theTrackMonitorTree = nullptr;
0502 }
0503 if (theHitMonitorTree != nullptr) {
0504 theHitMonitorTree->Write();
0505 delete theHitMonitorTree;
0506 theHitMonitorTree = nullptr;
0507 }
0508 theTrackHitMonitorIORootFile->Close();
0509 }
0510 }
0511
0512 bool HIPAlignmentAlgorithm::processHit1D(const AlignableDetOrUnitPtr& alidet,
0513 const Alignable* ali,
0514 const HIPAlignableSpecificParameters* alispecifics,
0515 const TrajectoryStateOnSurface& tsos,
0516 const TransientTrackingRecHit* hit,
0517 double hitwt) {
0518 static const unsigned int hitDim = 1;
0519 if (hitwt == 0.)
0520 return false;
0521
0522
0523 LocalPoint alvec = tsos.localPosition();
0524 AlgebraicVector pos(hitDim);
0525 pos[0] = alvec.x();
0526
0527
0528 AlgebraicSymMatrix ipcovmat(hitDim);
0529 ipcovmat[0][0] = tsos.localError().positionError().xx();
0530
0531
0532 AlgebraicVector coor(hitDim);
0533 coor[0] = hit->localPosition().x();
0534
0535 AlgebraicSymMatrix covmat(hitDim);
0536 covmat[0][0] = hit->localPositionError().xx();
0537
0538
0539 covmat = covmat + ipcovmat;
0540
0541
0542 double xpull = 0.;
0543 if (covmat[0][0] != 0.)
0544 xpull = (pos[0] - coor[0]) / sqrt(fabs(covmat[0][0]));
0545
0546
0547 AlignmentParameters* params = ali->alignmentParameters();
0548 HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
0549 uservar->datatype = theDataGroup;
0550
0551 AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
0552
0553 int npar = derivs2D.num_row();
0554 AlgebraicMatrix derivs(npar, hitDim, 0);
0555
0556 for (int ipar = 0; ipar < npar; ipar++) {
0557 for (unsigned int idim = 0; idim < hitDim; idim++) {
0558 derivs[ipar][idim] = derivs2D[ipar][idim];
0559 }
0560 }
0561
0562
0563 int ierr;
0564 covmat.invert(ierr);
0565 if (ierr != 0) {
0566 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::processHit1D"
0567 << "Cov. matrix inversion failed!";
0568 return false;
0569 }
0570
0571 double maxHitPullThreshold =
0572 (!theApplyCutsPerComponent ? defaultAlignableSpecs.maxHitPull : alispecifics->maxHitPull);
0573 bool useThisHit = (maxHitPullThreshold < 0.);
0574 useThisHit |= (fabs(xpull) < maxHitPullThreshold);
0575 if (!useThisHit) {
0576 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::processHit2D"
0577 << "Hit pull (x) " << xpull << " fails the cut " << maxHitPullThreshold;
0578 return false;
0579 }
0580
0581 AlgebraicMatrix covtmp(covmat);
0582 AlgebraicMatrix jtvjtmp(derivs * covtmp * derivs.T());
0583 AlgebraicSymMatrix thisjtvj(npar);
0584 AlgebraicVector thisjtve(npar);
0585 thisjtvj.assign(jtvjtmp);
0586 thisjtve = derivs * covmat * (pos - coor);
0587
0588 AlgebraicVector hitresidual(hitDim);
0589 hitresidual[0] = (pos[0] - coor[0]);
0590
0591 AlgebraicMatrix hitresidualT;
0592 hitresidualT = hitresidual.T();
0593
0594 uservar->jtvj += hitwt * thisjtvj;
0595 uservar->jtve += hitwt * thisjtve;
0596 uservar->nhit++;
0597
0598
0599 float thischi2;
0600 thischi2 = hitwt * (hitresidualT * covmat * hitresidual)[0];
0601
0602 if (verbose && (thischi2 / static_cast<float>(uservar->nhit)) > 10.)
0603 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::processHit1D] Added to Chi2 the number " << thischi2
0604 << " having " << uservar->nhit << " ndof"
0605 << ", X-resid " << hitresidual[0] << ", Cov^-1 matr (covmat):"
0606 << " [0][0] = " << covmat[0][0];
0607
0608 uservar->alichi2 += thischi2;
0609 uservar->alindof += hitDim;
0610
0611 return true;
0612 }
0613
0614 bool HIPAlignmentAlgorithm::processHit2D(const AlignableDetOrUnitPtr& alidet,
0615 const Alignable* ali,
0616 const HIPAlignableSpecificParameters* alispecifics,
0617 const TrajectoryStateOnSurface& tsos,
0618 const TransientTrackingRecHit* hit,
0619 double hitwt) {
0620 static const unsigned int hitDim = 2;
0621 if (hitwt == 0.)
0622 return false;
0623
0624
0625 LocalPoint alvec = tsos.localPosition();
0626 AlgebraicVector pos(hitDim);
0627 pos[0] = alvec.x();
0628 pos[1] = alvec.y();
0629
0630
0631 AlgebraicSymMatrix ipcovmat(hitDim);
0632 ipcovmat[0][0] = tsos.localError().positionError().xx();
0633 ipcovmat[1][1] = tsos.localError().positionError().yy();
0634 ipcovmat[0][1] = tsos.localError().positionError().xy();
0635
0636
0637 AlgebraicVector coor(hitDim);
0638 coor[0] = hit->localPosition().x();
0639 coor[1] = hit->localPosition().y();
0640
0641 AlgebraicSymMatrix covmat(hitDim);
0642 covmat[0][0] = hit->localPositionError().xx();
0643 covmat[1][1] = hit->localPositionError().yy();
0644 covmat[0][1] = hit->localPositionError().xy();
0645
0646
0647 covmat = covmat + ipcovmat;
0648
0649
0650 double xpull = 0.;
0651 double ypull = 0.;
0652 if (covmat[0][0] != 0.)
0653 xpull = (pos[0] - coor[0]) / sqrt(fabs(covmat[0][0]));
0654 if (covmat[1][1] != 0.)
0655 ypull = (pos[1] - coor[1]) / sqrt(fabs(covmat[1][1]));
0656
0657
0658 AlignmentParameters* params = ali->alignmentParameters();
0659 HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
0660 uservar->datatype = theDataGroup;
0661
0662 AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
0663
0664 int npar = derivs2D.num_row();
0665 AlgebraicMatrix derivs(npar, hitDim, 0);
0666
0667 for (int ipar = 0; ipar < npar; ipar++) {
0668 for (unsigned int idim = 0; idim < hitDim; idim++) {
0669 derivs[ipar][idim] = derivs2D[ipar][idim];
0670 }
0671 }
0672
0673
0674 int ierr;
0675 covmat.invert(ierr);
0676 if (ierr != 0) {
0677 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::processHit2D"
0678 << "Cov. matrix inversion failed!";
0679 return false;
0680 }
0681
0682 double maxHitPullThreshold =
0683 (!theApplyCutsPerComponent ? defaultAlignableSpecs.maxHitPull : alispecifics->maxHitPull);
0684 bool useThisHit = (maxHitPullThreshold < 0.);
0685 useThisHit |= (fabs(xpull) < maxHitPullThreshold && fabs(ypull) < maxHitPullThreshold);
0686 if (!useThisHit) {
0687 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::processHit2D"
0688 << "Hit pull (x,y) " << xpull << " , " << ypull << " fails the cut "
0689 << maxHitPullThreshold;
0690 return false;
0691 }
0692
0693 AlgebraicMatrix covtmp(covmat);
0694 AlgebraicMatrix jtvjtmp(derivs * covtmp * derivs.T());
0695 AlgebraicSymMatrix thisjtvj(npar);
0696 AlgebraicVector thisjtve(npar);
0697 thisjtvj.assign(jtvjtmp);
0698 thisjtve = derivs * covmat * (pos - coor);
0699
0700 AlgebraicVector hitresidual(hitDim);
0701 hitresidual[0] = (pos[0] - coor[0]);
0702 hitresidual[1] = (pos[1] - coor[1]);
0703
0704 AlgebraicMatrix hitresidualT;
0705 hitresidualT = hitresidual.T();
0706
0707 uservar->jtvj += hitwt * thisjtvj;
0708 uservar->jtve += hitwt * thisjtve;
0709 uservar->nhit++;
0710
0711
0712 float thischi2;
0713 thischi2 = hitwt * (hitresidualT * covmat * hitresidual)[0];
0714
0715 if (verbose && (thischi2 / static_cast<float>(uservar->nhit)) > 10.)
0716 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::processHit2D] Added to Chi2 the number " << thischi2
0717 << " having " << uservar->nhit << " ndof"
0718 << ", X-resid " << hitresidual[0] << ", Y-resid " << hitresidual[1]
0719 << ", Cov^-1 matr (covmat):"
0720 << " [0][0] = " << covmat[0][0] << " [0][1] = " << covmat[0][1]
0721 << " [1][0] = " << covmat[1][0] << " [1][1] = " << covmat[1][1];
0722
0723 uservar->alichi2 += thischi2;
0724 uservar->alindof += hitDim;
0725
0726 return true;
0727 }
0728
0729
0730 void HIPAlignmentAlgorithm::run(const edm::EventSetup& setup, const EventInfo& eventInfo) {
0731 if (isCollector)
0732 return;
0733
0734 TrajectoryStateCombiner tsoscomb;
0735
0736 m_datatype = theDataGroup;
0737
0738
0739 const ConstTrajTrackPairCollection& tracks = eventInfo.trajTrackPairs();
0740 for (ConstTrajTrackPairCollection::const_iterator it = tracks.begin(); it != tracks.end(); ++it) {
0741 const Trajectory* traj = (*it).first;
0742 const reco::Track* track = (*it).second;
0743
0744 float pt = track->pt();
0745 float eta = track->eta();
0746 float phi = track->phi();
0747 float p = track->p();
0748 float chi2n = track->normalizedChi2();
0749 int nhit = track->numberOfValidHits();
0750 float d0 = track->d0();
0751 float dz = track->dz();
0752
0753 int nhpxb = track->hitPattern().numberOfValidPixelBarrelHits();
0754 int nhpxf = track->hitPattern().numberOfValidPixelEndcapHits();
0755 int nhtib = track->hitPattern().numberOfValidStripTIBHits();
0756 int nhtob = track->hitPattern().numberOfValidStripTOBHits();
0757 int nhtid = track->hitPattern().numberOfValidStripTIDHits();
0758 int nhtec = track->hitPattern().numberOfValidStripTECHits();
0759
0760 if (verbose)
0761 edm::LogInfo("Alignment") << "New track pt,eta,phi,chi2n,hits: " << pt << "," << eta << "," << phi << "," << chi2n
0762 << "," << nhit;
0763
0764 double ihitwt = 1;
0765 double trkwt = 1;
0766 if (trackWt) {
0767 trkwt = Scale;
0768
0769 if (uniEta)
0770 trkwt *= theEtaFormula->Eval(fabs(eta));
0771 }
0772 if (trackPs) {
0773 double r = gRandom->Rndm();
0774 if (trkwt < r)
0775 continue;
0776 } else if (trackWt)
0777 ihitwt = trkwt;
0778
0779
0780 {
0781 theMonitorConfig.trackmonitorvars.m_Nhits.push_back(nhit);
0782 theMonitorConfig.trackmonitorvars.m_Pt.push_back(pt);
0783 theMonitorConfig.trackmonitorvars.m_P.push_back(p);
0784 theMonitorConfig.trackmonitorvars.m_Eta.push_back(eta);
0785 theMonitorConfig.trackmonitorvars.m_Phi.push_back(phi);
0786 theMonitorConfig.trackmonitorvars.m_Chi2n.push_back(chi2n);
0787 theMonitorConfig.trackmonitorvars.m_nhPXB.push_back(nhpxb);
0788 theMonitorConfig.trackmonitorvars.m_nhPXF.push_back(nhpxf);
0789 theMonitorConfig.trackmonitorvars.m_nhTIB.push_back(nhtib);
0790 theMonitorConfig.trackmonitorvars.m_nhTOB.push_back(nhtob);
0791 theMonitorConfig.trackmonitorvars.m_nhTID.push_back(nhtid);
0792 theMonitorConfig.trackmonitorvars.m_nhTEC.push_back(nhtec);
0793 theMonitorConfig.trackmonitorvars.m_d0.push_back(d0);
0794 theMonitorConfig.trackmonitorvars.m_dz.push_back(dz);
0795 theMonitorConfig.trackmonitorvars.m_wt.push_back(ihitwt);
0796 }
0797
0798 std::vector<const TransientTrackingRecHit*> hitvec;
0799 std::vector<TrajectoryStateOnSurface> tsosvec;
0800
0801
0802 std::vector<TrajectoryMeasurement> measurements = traj->measurements();
0803 for (std::vector<TrajectoryMeasurement>::iterator im = measurements.begin(); im != measurements.end(); ++im) {
0804 TrajectoryMeasurement meas = *im;
0805
0806
0807
0808 const TransientTrackingRecHit* hit = &(*meas.recHit());
0809
0810 if (hit->isValid() && theAlignableDetAccessor->detAndSubdetInMap(hit->geographicalId())) {
0811
0812
0813
0814
0815
0816
0817 if (eventInfo.clusterValueMap()) {
0818
0819
0820
0821 AlignmentClusterFlag myflag;
0822
0823 int subDet = hit->geographicalId().subdetId();
0824
0825 const TrackingRecHit* rechit = hit->hit();
0826 if (subDet > 2) {
0827 const std::type_info& type = typeid(*rechit);
0828
0829 if (type == typeid(SiStripRecHit1D)) {
0830 const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(rechit);
0831 if (stripHit1D) {
0832 SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
0833
0834 myflag = (*eventInfo.clusterValueMap())[stripclust];
0835 } else
0836 edm::LogError("HIPAlignmentAlgorithm")
0837 << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
0838 << "TypeId of the RecHit: " << className(*hit);
0839 }
0840 else if (type == typeid(SiStripRecHit2D)) {
0841 const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(rechit);
0842 if (stripHit2D) {
0843 SiStripRecHit2D::ClusterRef stripclust(stripHit2D->cluster());
0844
0845 myflag = (*eventInfo.clusterValueMap())[stripclust];
0846 } else
0847 edm::LogError("HIPAlignmentAlgorithm")
0848 << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
0849 << "TypeId of the TTRH: " << className(*hit);
0850 }
0851 }
0852 else {
0853 const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(rechit);
0854 if (pixelhit) {
0855 SiPixelClusterRefNew pixelclust(pixelhit->cluster());
0856
0857 myflag = (*eventInfo.clusterValueMap())[pixelclust];
0858 } else
0859 edm::LogError("HIPAlignmentAlgorithm")
0860 << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Pixel RecHit failed! "
0861 << "TypeId of the TTRH: " << className(*hit);
0862 }
0863 if (!myflag.isTaken())
0864 continue;
0865 }
0866
0867
0868 TrajectoryStateOnSurface tsos = tsoscomb.combine(meas.forwardPredictedState(), meas.backwardPredictedState());
0869
0870 if (tsos.isValid()) {
0871 hitvec.push_back(hit);
0872 tsosvec.push_back(tsos);
0873 }
0874
0875 }
0876 }
0877
0878
0879 std::vector<AlignableDetOrUnitPtr> alidetvec = theAlignableDetAccessor->alignablesFromHits(hitvec);
0880
0881
0882 CompositeAlignmentParameters aap = theAlignmentParameterStore->selectParameters(alidetvec);
0883
0884 std::vector<TrajectoryStateOnSurface>::const_iterator itsos = tsosvec.begin();
0885 std::vector<const TransientTrackingRecHit*>::const_iterator ihit = hitvec.begin();
0886
0887
0888 while (itsos != tsosvec.end()) {
0889
0890 const GeomDet* det = (*ihit)->det();
0891
0892 uint32_t nhitDim = (*ihit)->dimension();
0893
0894 AlignableDetOrUnitPtr alidet = theAlignableDetAccessor->alignableFromGeomDet(det);
0895
0896
0897 Alignable* ali = aap.alignableFromAlignableDet(alidet);
0898
0899 if (ali != nullptr) {
0900 const HIPAlignableSpecificParameters* alispecifics = findAlignableSpecs(ali);
0901 const TrajectoryStateOnSurface& tsos = *itsos;
0902
0903
0904
0905
0906
0907 double mom_x = tsos.localDirection().x();
0908 double mom_y = tsos.localDirection().y();
0909 double mom_z = tsos.localDirection().z();
0910 double sin_theta = TMath::Abs(mom_z) / sqrt(pow(mom_x, 2) + pow(mom_y, 2) + pow(mom_z, 2));
0911 double angle = TMath::ASin(sin_theta);
0912 double alihitwt = ihitwt;
0913
0914
0915 if (IsCollision) {
0916 if (angle > col_cut)
0917 alihitwt = 0;
0918 } else {
0919 if (angle < cos_cut)
0920 alihitwt = 0;
0921 }
0922
0923
0924 theMonitorConfig.hitmonitorvars.m_angle = angle;
0925 theMonitorConfig.hitmonitorvars.m_sinTheta = sin_theta;
0926 theMonitorConfig.hitmonitorvars.m_detId = ali->id();
0927
0928
0929 if ((*ihit)->hit() != nullptr) {
0930 const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>((*ihit)->hit());
0931 if (pixhit != nullptr) {
0932 theMonitorConfig.hitmonitorvars.m_hasHitProb = pixhit->hasFilledProb();
0933 if (theMonitorConfig.hitmonitorvars.m_hasHitProb) {
0934
0935 theMonitorConfig.hitmonitorvars.m_probXY = pixhit->probabilityXY();
0936 theMonitorConfig.hitmonitorvars.m_probQ = pixhit->probabilityQ();
0937 theMonitorConfig.hitmonitorvars.m_rawQualityWord = pixhit->rawQualityWord();
0938 if (alispecifics->applyPixelProbCut) {
0939 bool probXYgood = (theMonitorConfig.hitmonitorvars.m_probXY >= alispecifics->minPixelProbXY &&
0940 theMonitorConfig.hitmonitorvars.m_probXY <= alispecifics->maxPixelProbXY);
0941 bool probQgood = (theMonitorConfig.hitmonitorvars.m_probQ >= alispecifics->minPixelProbQ &&
0942 theMonitorConfig.hitmonitorvars.m_probQ <= alispecifics->maxPixelProbQ);
0943 bool probXYQgood;
0944 if (alispecifics->usePixelProbXYOrProbQ)
0945 probXYQgood = (probXYgood || probQgood);
0946 else
0947 probXYQgood = (probXYgood && probQgood);
0948 if (!probXYQgood)
0949 alihitwt = 0;
0950 }
0951 }
0952 }
0953 }
0954
0955 theMonitorConfig.hitmonitorvars.m_hitwt = alihitwt;
0956 bool hitProcessed = false;
0957 switch (nhitDim) {
0958 case 1:
0959 hitProcessed = processHit1D(alidet, ali, alispecifics, tsos, *ihit, alihitwt);
0960 break;
0961 case 2:
0962 hitProcessed = processHit2D(alidet, ali, alispecifics, tsos, *ihit, alihitwt);
0963 break;
0964 default:
0965 edm::LogError("HIPAlignmentAlgorithm")
0966 << "ERROR in <HIPAlignmentAlgorithm::run>: Number of hit dimensions = " << nhitDim
0967 << " is not supported!" << std::endl;
0968 break;
0969 }
0970 if (theMonitorConfig.fillTrackHitMonitoring && theMonitorConfig.checkNhits() && hitProcessed)
0971 theMonitorConfig.hitmonitorvars.fill();
0972 }
0973
0974 itsos++;
0975 ihit++;
0976 }
0977 }
0978
0979
0980 if (theMonitorConfig.fillTrackMonitoring && theMonitorConfig.checkNevents())
0981 theMonitorConfig.trackmonitorvars.fill();
0982 }
0983
0984
0985 int HIPAlignmentAlgorithm::readIterationFile(std::string filename) {
0986 int result;
0987
0988 std::ifstream inIterFile(filename.c_str(), std::ios::in);
0989 if (!inIterFile) {
0990 edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] ERROR! "
0991 << "Unable to open Iteration file";
0992 result = -1;
0993 } else {
0994 inIterFile >> result;
0995 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] "
0996 << "Read last iteration number from file: " << result;
0997 }
0998 inIterFile.close();
0999
1000 return result;
1001 }
1002
1003
1004 void HIPAlignmentAlgorithm::writeIterationFile(std::string filename, int iter) {
1005 std::ofstream outIterFile((filename.c_str()), std::ios::out);
1006 if (!outIterFile)
1007 edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::writeIterationFile] ERROR: Unable to write Iteration file";
1008 else {
1009 outIterFile << iter;
1010 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::writeIterationFile] writing iteration number to file: "
1011 << iter;
1012 }
1013 outIterFile.close();
1014 }
1015
1016
1017
1018 void HIPAlignmentAlgorithm::setAlignmentPositionError(void) {
1019
1020 if (!theApplyAPE) {
1021 edm::LogInfo("Alignment") << "[HIPAlignmentAlgorithm::setAlignmentPositionError] No APE applied";
1022 return;
1023 }
1024
1025 edm::LogInfo("Alignment") << "[HIPAlignmentAlgorithm::setAlignmentPositionError] Apply APE!";
1026
1027 double apeSPar[3], apeRPar[3];
1028 for (const auto& alipars : theAPEParameters) {
1029 const auto& alignables = alipars.first;
1030 const auto& pars = alipars.second;
1031
1032 apeSPar[0] = pars[0];
1033 apeSPar[1] = pars[1];
1034 apeSPar[2] = pars[2];
1035 apeRPar[0] = pars[3];
1036 apeRPar[1] = pars[4];
1037 apeRPar[2] = pars[5];
1038
1039 int function = pars[6];
1040
1041
1042 printf("APE: %u alignables\n", (unsigned int)alignables.size());
1043 for (int i = 0; i < 21; ++i) {
1044 double apelinstest = calcAPE(apeSPar, i, 0);
1045 double apeexpstest = calcAPE(apeSPar, i, 1);
1046 double apestepstest = calcAPE(apeSPar, i, 2);
1047 double apelinrtest = calcAPE(apeRPar, i, 0);
1048 double apeexprtest = calcAPE(apeRPar, i, 1);
1049 double apesteprtest = calcAPE(apeRPar, i, 2);
1050 printf("APE: iter slin sexp sstep rlin rexp rstep: %5d %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f \n",
1051 i,
1052 apelinstest,
1053 apeexpstest,
1054 apestepstest,
1055 apelinrtest,
1056 apeexprtest,
1057 apesteprtest);
1058 }
1059
1060
1061 double apeshift = calcAPE(apeSPar, theIteration, function);
1062 double aperot = calcAPE(apeRPar, theIteration, function);
1063 theAlignmentParameterStore->setAlignmentPositionError(alignables, apeshift, aperot);
1064 }
1065 }
1066
1067
1068
1069 double HIPAlignmentAlgorithm::calcAPE(double* par, int iter, int function) {
1070 double diter = (double)iter;
1071 if (function == 0)
1072 return std::max(par[1], par[0] + ((par[1] - par[0]) / par[2]) * diter);
1073 else if (function == 1)
1074 return std::max(0., par[0] * (exp(-pow(diter, par[1]) / par[2])));
1075 else if (function == 2) {
1076 int ipar2 = (int)par[2];
1077 int step = iter / ipar2;
1078 double dstep = (double)step;
1079 return std::max(0., par[0] - par[1] * dstep);
1080 } else
1081 assert(false);
1082 }
1083
1084
1085
1086 void HIPAlignmentAlgorithm::bookRoot(void) {
1087
1088 if (doTrackHitMonitoring && !isCollector) {
1089 theTrackHitMonitorIORootFile = TFile::Open(theMonitorConfig.outfile.c_str(), "update");
1090 theTrackHitMonitorIORootFile->cd();
1091
1092 if (theMonitorConfig.fillTrackMonitoring) {
1093 TString tname = Form("T1_%i", theIteration);
1094 theTrackMonitorTree = new TTree(tname, "Eventwise tree");
1095
1096
1097 theTrackMonitorTree->Branch("DataType", &m_datatype);
1098 theMonitorConfig.trackmonitorvars.setTree(theTrackMonitorTree);
1099 theMonitorConfig.trackmonitorvars.bookBranches();
1100 }
1101
1102 if (theMonitorConfig.fillTrackHitMonitoring) {
1103 TString tname_hit = Form("T1_hit_%i", theIteration);
1104 theHitMonitorTree = new TTree(tname_hit, "Hitwise tree");
1105 theHitMonitorTree->Branch("DataType", &m_datatype);
1106 theMonitorConfig.hitmonitorvars.setTree(theHitMonitorTree);
1107 theMonitorConfig.hitmonitorvars.bookBranches();
1108 }
1109 }
1110
1111
1112 if (isCollector) {
1113 TString tname = Form("T2_%i", theIteration);
1114 theAlignablesMonitorIORootFile = TFile::Open(outfile2.c_str(), "update");
1115 theAlignablesMonitorIORootFile->cd();
1116 theAlignablesMonitorTree = new TTree(tname, "Alignablewise tree");
1117 theAlignablesMonitorTree->Branch("Id", &m2_Id, "Id/i");
1118 theAlignablesMonitorTree->Branch("ObjId", &m2_ObjId, "ObjId/I");
1119 theAlignablesMonitorTree->Branch("Nhit", &m2_Nhit);
1120 theAlignablesMonitorTree->Branch("DataType", &m2_datatype);
1121 theAlignablesMonitorTree->Branch("Type", &m2_Type);
1122 theAlignablesMonitorTree->Branch("Layer", &m2_Layer);
1123 theAlignablesMonitorTree->Branch("Xpos", &m2_Xpos);
1124 theAlignablesMonitorTree->Branch("Ypos", &m2_Ypos);
1125 theAlignablesMonitorTree->Branch("Zpos", &m2_Zpos);
1126 theAlignablesMonitorTree->Branch("DeformationsType", &m2_dtype, "DeformationsType/I");
1127 theAlignablesMonitorTree->Branch("NumDeformations", &m2_nsurfdef);
1128 theAlignablesMonitorTree->Branch("Deformations", &m2_surfDef);
1129 }
1130
1131
1132 if (!theLevels.empty()) {
1133 TString tname = Form("T3_%i", theIteration);
1134 theSurveyIORootFile = TFile::Open(ssurveyfile.c_str(), "update");
1135 theSurveyIORootFile->cd();
1136 theSurveyTree = new TTree(tname, "Survey Tree");
1137 theSurveyTree->Branch("Id", &m3_Id, "Id/i");
1138 theSurveyTree->Branch("ObjId", &m3_ObjId, "ObjId/I");
1139 theSurveyTree->Branch("Par", &m3_par, "Par[6]/F");
1140 edm::LogInfo("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Survey trees booked.";
1141 }
1142 edm::LogInfo("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Root trees booked.";
1143 }
1144
1145
1146
1147 void HIPAlignmentAlgorithm::fillAlignablesMonitor(const edm::EventSetup& iSetup) {
1148 if (theAlignablesMonitorIORootFile == (TFile*)nullptr)
1149 return;
1150 using std::setw;
1151 theAlignablesMonitorIORootFile->cd();
1152
1153 int naligned = 0;
1154
1155
1156 const TrackerTopology* const tTopo = &iSetup.getData(topoToken2_);
1157
1158 for (const auto& ali : theAlignables) {
1159 AlignmentParameters* dap = ali->alignmentParameters();
1160
1161
1162 if (dap->isValid()) {
1163
1164 HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(dap->userVariables());
1165 m2_Nhit = uservar->nhit;
1166 m2_datatype = uservar->datatype;
1167
1168
1169 std::pair<int, int> tl = theAlignmentParameterStore->typeAndLayer(ali, tTopo);
1170 m2_Type = tl.first;
1171 m2_Layer = tl.second;
1172
1173
1174 m2_Id = ali->id();
1175 m2_ObjId = ali->alignableObjectId();
1176
1177
1178 GlobalPoint pos = ali->surface().position();
1179 m2_Xpos = pos.x();
1180 m2_Ypos = pos.y();
1181 m2_Zpos = pos.z();
1182
1183 m2_surfDef.clear();
1184 {
1185 std::vector<std::pair<int, SurfaceDeformation*>> dali_id_pairs;
1186 SurfaceDeformation* dali_obj = nullptr;
1187 SurfaceDeformationFactory::Type dtype = SurfaceDeformationFactory::kNoDeformations;
1188 std::vector<double> dali;
1189 if (1 == ali->surfaceDeformationIdPairs(dali_id_pairs)) {
1190 dali_obj = dali_id_pairs[0].second;
1191 dali = dali_obj->parameters();
1192 dtype = (SurfaceDeformationFactory::Type)dali_obj->type();
1193 }
1194 for (auto& dit : dali)
1195 m2_surfDef.push_back((float)dit);
1196 m2_dtype = dtype;
1197 m2_nsurfdef = m2_surfDef.size();
1198 }
1199
1200 if (verbose) {
1201 AlgebraicVector pars = dap->parameters();
1202 edm::LogVerbatim("Alignment") << "------------------------------------------------------------------------\n"
1203 << " ALIGNABLE: " << setw(6) << naligned << '\n'
1204 << "hits: " << setw(4) << m2_Nhit << " type: " << setw(4) << m2_Type
1205 << " layer: " << setw(4) << m2_Layer << " id: " << setw(4) << m2_Id
1206 << " objId: " << setw(4) << m2_ObjId << '\n'
1207 << std::fixed << std::setprecision(5) << "x,y,z: " << setw(12) << m2_Xpos
1208 << setw(12) << m2_Ypos << setw(12) << m2_Zpos
1209 << "\nDeformations type, nDeformations: " << setw(12) << m2_dtype << setw(12)
1210 << m2_nsurfdef << '\n'
1211 << "params: " << setw(12) << pars[0] << setw(12) << pars[1] << setw(12) << pars[2]
1212 << setw(12) << pars[3] << setw(12) << pars[4] << setw(12) << pars[5];
1213 }
1214
1215 naligned++;
1216 if (theAlignablesMonitorTree != nullptr)
1217 theAlignablesMonitorTree->Fill();
1218 }
1219 }
1220 }
1221
1222
1223 bool HIPAlignmentAlgorithm::calcParameters(Alignable* ali, int setDet, double start, double step) {
1224 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::calcParameters"
1225 << "Begin: Processing detector " << ali->id();
1226
1227
1228 AlignmentParameters* par = ali->alignmentParameters();
1229 const HIPAlignableSpecificParameters* alispecifics = findAlignableSpecs(ali);
1230
1231 HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(par->userVariables());
1232 int nhit = uservar->nhit;
1233
1234
1235
1236
1237
1238 int minHitThreshold = (!theApplyCutsPerComponent ? defaultAlignableSpecs.minNHits : alispecifics->minNHits);
1239 if (!isCollector)
1240 minHitThreshold = 1;
1241 if (setDet == 0 && nhit < minHitThreshold) {
1242 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::calcParameters"
1243 << "Skipping detector " << ali->id() << " because number of hits = " << nhit
1244 << " <= " << minHitThreshold;
1245 par->setValid(false);
1246 return false;
1247 }
1248
1249 AlgebraicSymMatrix jtvj = uservar->jtvj;
1250 AlgebraicVector jtve = uservar->jtve;
1251
1252
1253 int npar = jtve.num_row();
1254 AlgebraicVector params(npar);
1255 AlgebraicVector paramerr(npar);
1256 AlgebraicSymMatrix cov(npar);
1257
1258
1259 if (isCollector) {
1260 if (setDet == 0) {
1261 int ierr;
1262 AlgebraicSymMatrix jtvjinv = jtvj.inverse(ierr);
1263 if (ierr != 0) {
1264 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::calcParameters"
1265 << "Matrix inversion failed!";
1266 return false;
1267 }
1268 params = -(jtvjinv * jtve);
1269 cov = jtvjinv;
1270
1271 double minRelErrThreshold =
1272 (!theApplyCutsPerComponent ? defaultAlignableSpecs.minRelParError : alispecifics->minRelParError);
1273 double maxRelErrThreshold =
1274 (!theApplyCutsPerComponent ? defaultAlignableSpecs.maxRelParError : alispecifics->maxRelParError);
1275 for (int i = 0; i < npar; i++) {
1276 double relerr = 0;
1277 if (fabs(cov[i][i]) > 0.)
1278 paramerr[i] = sqrt(fabs(cov[i][i]));
1279 else
1280 paramerr[i] = params[i];
1281 if (params[i] != 0.)
1282 relerr = fabs(paramerr[i] / params[i]);
1283 if ((maxRelErrThreshold >= 0. && relerr > maxRelErrThreshold) || relerr < minRelErrThreshold) {
1284 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::calcParameters"
1285 << "RelError = " << relerr << " > " << maxRelErrThreshold << " or < "
1286 << minRelErrThreshold << ". Setting param = paramerr = 0 for component " << i;
1287 params[i] = 0;
1288 paramerr[i] = 0;
1289 }
1290 }
1291 } else {
1292 if (params.num_row() != 1) {
1293 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::calcParameters"
1294 << "For scanning, please only turn on one parameter! check common_cff_py.txt";
1295 return false;
1296 }
1297 if (theIteration == 1)
1298 params[0] = start;
1299 else
1300 params[0] = step;
1301 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::calcParameters"
1302 << "Parameters = " << params;
1303 }
1304 }
1305
1306 uservar->alipar = params;
1307 uservar->alierr = paramerr;
1308
1309 AlignmentParameters* parnew = par->cloneFromSelected(params, cov);
1310 ali->setAlignmentParameters(parnew);
1311 parnew->setValid(true);
1312
1313 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::calcParameters"
1314 << "End: Processing detector " << ali->id();
1315
1316 return true;
1317 }
1318
1319
1320 void HIPAlignmentAlgorithm::collector(void) {
1321 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1322 << "Called for iteration " << theIteration;
1323
1324 std::vector<std::string> monitorFileList;
1325 HIPUserVariablesIORoot HIPIO;
1326
1327 typedef int pawt_t;
1328 std::unordered_map<Alignable*, std::unordered_map<int, pawt_t>> ali_datatypecountpair_map;
1329 if (rewgtPerAli) {
1330 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1331 << "Per-alignable reweighting is turned on. Iterating over the parallel jobs to sum "
1332 "number of hits per alignable for each data type.";
1333
1334 for (int ijob = 1; ijob <= theCollectorNJobs; ijob++) {
1335 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1336 << "Pre-reading uservar for job " << ijob;
1337
1338 std::string str = std::to_string(ijob);
1339 std::string uvfile = theCollectorPath + "/job" + str + "/" + suvarfilecore;
1340
1341 std::vector<AlignmentUserVariables*> uvarvec =
1342 HIPIO.readHIPUserVariables(theAlignables, uvfile.c_str(), theIteration, ioerr);
1343 if (uvarvec.size() != theAlignables.size())
1344 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1345 << "Number of alignables = " << theAlignables.size()
1346 << " is not the same as number of user variables = " << uvarvec.size()
1347 << ". A mismatch might occur!";
1348 if (ioerr != 0) {
1349 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1350 << "Could not read user variable files for job " << ijob << " in iteration "
1351 << theIteration;
1352 continue;
1353 }
1354 std::vector<AlignmentUserVariables*>::const_iterator iuvar =
1355 uvarvec.begin();
1356 for (const auto& ali : theAlignables) {
1357
1358
1359 HIPUserVariables* uvar = dynamic_cast<HIPUserVariables*>(*iuvar);
1360 if (uvar != nullptr) {
1361 int alijobdtype = uvar->datatype;
1362 pawt_t alijobnhits = uvar->nhit;
1363 if (ali_datatypecountpair_map.find(ali) == ali_datatypecountpair_map.end()) {
1364 std::unordered_map<int, pawt_t> newmap;
1365 ali_datatypecountpair_map[ali] = newmap;
1366 ali_datatypecountpair_map[ali][alijobdtype] = alijobnhits;
1367 } else {
1368 std::unordered_map<int, pawt_t>& theMap = ali_datatypecountpair_map[ali];
1369 if (theMap.find(alijobdtype) == theMap.end())
1370 theMap[alijobdtype] = alijobnhits;
1371 else
1372 theMap[alijobdtype] += alijobnhits;
1373 }
1374 delete uvar;
1375 }
1376 iuvar++;
1377 }
1378 }
1379 }
1380
1381 for (int ijob = 1; ijob <= theCollectorNJobs; ijob++) {
1382 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1383 << "Reading uservar for job " << ijob;
1384
1385 std::string str = std::to_string(ijob);
1386 std::string uvfile = theCollectorPath + "/job" + str + "/" + suvarfilecore;
1387
1388 std::vector<AlignmentUserVariables*> uvarvec =
1389 HIPIO.readHIPUserVariables(theAlignables, uvfile.c_str(), theIteration, ioerr);
1390 if (uvarvec.size() != theAlignables.size())
1391 edm::LogWarning("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1392 << "Number of alignables = " << theAlignables.size()
1393 << " is not the same as number of user variables = " << uvarvec.size()
1394 << ". A mismatch might occur!";
1395
1396 if (ioerr != 0) {
1397 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1398 << "Could not read user variable files for job " << ijob << " in iteration "
1399 << theIteration;
1400 continue;
1401 }
1402
1403
1404 std::vector<AlignmentUserVariables*> uvarvecadd;
1405 std::vector<AlignmentUserVariables*>::const_iterator iuvarnew = uvarvec.begin();
1406 for (const auto& ali : theAlignables) {
1407 AlignmentParameters* ap = ali->alignmentParameters();
1408
1409 HIPUserVariables* uvarold = dynamic_cast<HIPUserVariables*>(ap->userVariables());
1410 HIPUserVariables* uvarnew = dynamic_cast<HIPUserVariables*>(*iuvarnew);
1411
1412 HIPUserVariables* uvar = uvarold->clone();
1413 uvar->datatype =
1414 theDataGroup;
1415
1416 if (uvarnew != nullptr) {
1417 double peraliwgt = 1;
1418 if (rewgtPerAli) {
1419 int alijobdtype = uvarnew->datatype;
1420 if (ali_datatypecountpair_map.find(ali) != ali_datatypecountpair_map.end() &&
1421 ali_datatypecountpair_map[ali].find(alijobdtype) != ali_datatypecountpair_map[ali].end()) {
1422 peraliwgt = ali_datatypecountpair_map[ali][alijobdtype];
1423 unsigned int nNonZeroTypes = 0;
1424 pawt_t sumwgts = 0;
1425 for (auto it = ali_datatypecountpair_map[ali].cbegin(); it != ali_datatypecountpair_map[ali].cend(); ++it) {
1426 sumwgts += it->second;
1427 if (it->second != pawt_t(0))
1428 nNonZeroTypes++;
1429 }
1430 edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1431 << "Reweighting detector " << ali->id() << " / " << ali->alignableObjectId()
1432 << " for data type " << alijobdtype << " by " << sumwgts << "/" << peraliwgt
1433 << "/" << nNonZeroTypes;
1434 peraliwgt = ((nNonZeroTypes == 0 || peraliwgt == double(0))
1435 ? double(1)
1436 : double((double(sumwgts)) / peraliwgt / (double(nNonZeroTypes))));
1437 } else if (ali_datatypecountpair_map.find(ali) != ali_datatypecountpair_map.end())
1438 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1439 << "Could not find data type " << alijobdtype << " for detector " << ali->id()
1440 << " / " << ali->alignableObjectId();
1441 else
1442 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::collector"
1443 << "Could not find detector " << ali->id() << " / " << ali->alignableObjectId()
1444 << " in the map ali_datatypecountpair_map";
1445 }
1446
1447 uvar->nhit = (uvarold->nhit) + (uvarnew->nhit);
1448 uvar->jtvj = (uvarold->jtvj) + peraliwgt * (uvarnew->jtvj);
1449 uvar->jtve = (uvarold->jtve) + peraliwgt * (uvarnew->jtve);
1450 uvar->alichi2 = (uvarold->alichi2) + peraliwgt * (uvarnew->alichi2);
1451 uvar->alindof = (uvarold->alindof) + (uvarnew->alindof);
1452
1453 delete uvarnew;
1454 }
1455
1456 uvarvecadd.push_back(uvar);
1457 iuvarnew++;
1458 }
1459
1460 theAlignmentParameterStore->attachUserVariables(theAlignables, uvarvecadd, ioerr);
1461
1462
1463 if (doTrackHitMonitoring) {
1464 uvfile = theCollectorPath + "/job" + str + "/" + theMonitorConfig.outfilecore;
1465 monitorFileList.push_back(uvfile);
1466 }
1467 }
1468
1469
1470 if (doTrackHitMonitoring)
1471 collectMonitorTrees(monitorFileList);
1472 }
1473
1474
1475 void HIPAlignmentAlgorithm::collectMonitorTrees(const std::vector<std::string>& filenames) {
1476 if (!doTrackHitMonitoring)
1477 return;
1478 if (!isCollector)
1479 throw cms::Exception("LogicError") << "[HIPAlignmentAlgorithm::collectMonitorTrees] Called in non-collector mode."
1480 << std::endl;
1481
1482 TString theTrackMonitorTreeName = Form("T1_%i", theIteration);
1483 TString theHitMonitorTreeName = Form("T1_hit_%i", theIteration);
1484
1485 std::vector<TFile*> finputlist;
1486 TList* eventtrees = new TList;
1487 TList* hittrees = new TList;
1488 for (std::string const& filename : filenames) {
1489 TFile* finput = TFile::Open(filename.c_str(), "read");
1490 if (finput != nullptr) {
1491 TTree* tmptree;
1492 if (theMonitorConfig.fillTrackMonitoring) {
1493 tmptree = nullptr;
1494 tmptree = (TTree*)finput->Get(theTrackMonitorTreeName);
1495 if (tmptree != nullptr)
1496 eventtrees->Add(tmptree);
1497 }
1498 if (theMonitorConfig.fillTrackHitMonitoring) {
1499 tmptree = nullptr;
1500 tmptree = (TTree*)finput->Get(theHitMonitorTreeName);
1501 if (tmptree != nullptr)
1502 hittrees->Add((TTree*)finput->Get(theHitMonitorTreeName));
1503 }
1504 finputlist.push_back(finput);
1505 }
1506 }
1507
1508 if (theTrackHitMonitorIORootFile != nullptr) {
1509 edm::LogError("Alignment") << "@SUB=HIPAlignmentAlgorithm::collectMonitorTrees"
1510 << "Monitor file is already open while it is not supposed to be!";
1511 delete theTrackMonitorTree;
1512 theTrackMonitorTree = nullptr;
1513 delete theHitMonitorTree;
1514 theHitMonitorTree = nullptr;
1515 theTrackHitMonitorIORootFile->Close();
1516 }
1517 theTrackHitMonitorIORootFile = TFile::Open(theMonitorConfig.outfile.c_str(), "update");
1518 theTrackHitMonitorIORootFile->cd();
1519 if (eventtrees->GetSize() > 0)
1520 theTrackMonitorTree = TTree::MergeTrees(eventtrees);
1521 if (hittrees->GetSize() > 0)
1522 theHitMonitorTree = TTree::MergeTrees(hittrees);
1523
1524
1525 delete hittrees;
1526 delete eventtrees;
1527 for (TFile*& finput : finputlist)
1528 finput->Close();
1529
1530
1531 if (theTrackMonitorTree != nullptr)
1532 theTrackMonitorTree->SetName(theTrackMonitorTreeName);
1533 if (theHitMonitorTree != nullptr)
1534 theHitMonitorTree->SetName(theHitMonitorTreeName);
1535 }
1536
1537
1538 HIPAlignableSpecificParameters* HIPAlignmentAlgorithm::findAlignableSpecs(const Alignable* ali) {
1539 if (ali != nullptr) {
1540 for (std::vector<HIPAlignableSpecificParameters>::iterator it = theAlignableSpecifics.begin();
1541 it != theAlignableSpecifics.end();
1542 it++) {
1543 if (it->matchAlignable(ali))
1544 return &(*it);
1545 }
1546 edm::LogInfo("Alignment") << "[HIPAlignmentAlgorithm::findAlignableSpecs] Alignment object with id " << ali->id()
1547 << " / " << ali->alignableObjectId() << " could not be found. Returning default.";
1548 }
1549 return &defaultAlignableSpecs;
1550 }