Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:23 

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 // Constructor ----------------------------------------------------------------
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   // parse parameters
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;  //Alignablewise tree
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   // parameters for APE
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   // for collector mode (parallel processing)
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 // Call at beginning of job ---------------------------------------------------
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           //sparameterfile.replace(sparameterfile.end()-5, sparameterfile.end(),iovapp);
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   // accessor Det->AlignableDet
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   // set alignmentParameterStore
0183   theAlignmentParameterStore = store;
0184 
0185   // get alignables
0186   theAlignables = theAlignmentParameterStore->alignables();
0187 
0188   // Config flags that specify different detectors
0189   {
0190     AlignmentParameterSelector selector(tracker, muon, extras);
0191 
0192     // APE parameters, clear if necessary
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);  // c.f. note in calcAPE
0222         else if (function == std::string("exponential"))
0223           apeSPar.push_back(1);  // c.f. note in calcAPE
0224         else if (function == std::string("step"))
0225           apeSPar.push_back(2);  // c.f. note in calcAPE
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     // Relative error per component instead of overall relative error
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 // Call at new loop -------------------------------------------------------------
0297 void HIPAlignmentAlgorithm::startNewLoop(void) {
0298   edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0299                             << "Begin";
0300 
0301   // iterate over all alignables and attach user variables
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   // try to read in alignment parameters from a previous iteration
0310   AlignablePositions theAlignablePositionsFromFile =
0311       theIO.readAlignableAbsolutePositions(theAlignables, salignedfile.c_str(), -1, ioerr);
0312   int numAlignablesFromFile = theAlignablePositionsFromFile.size();
0313   if (numAlignablesFromFile == 0) {  // file not there: first iteration
0314     // set iteration number to 1 when needed
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 {  // there have been previous iterations
0322     edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0323                               << "Alignables Read " << numAlignablesFromFile;
0324 
0325     // get iteration number from file
0326     theIteration = readIterationFile(siterationfile);
0327     // Where is the target for this?
0328     theIO.readAlignableAbsolutePositions(theAlignables, salignedfile.c_str(), theIteration, ioerr);
0329 
0330     // increase iteration
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     // now apply psotions of file from prev iteration
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   // book root trees
0347   bookRoot();
0348 
0349   // set alignment position error
0350   setAlignmentPositionError();
0351 
0352   // run collector job if we are in parallel mode
0353   if (isCollector)
0354     collector();
0355 
0356   edm::LogInfo("Alignment") << "@SUB=HIPAlignmentAlgorithm::startNewLoop"
0357                             << "End";
0358 }
0359 
0360 // Call at end of job ---------------------------------------------------------
0361 void HIPAlignmentAlgorithm::terminate(const edm::EventSetup& iSetup) {
0362   edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Terminating";
0363 
0364   // calculating survey residuals
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       // get position
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)) {  // FIXME: Needs revision for hardcoded consts
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             // variable for tree
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   // write user variables
0413   HIPUserVariablesIORoot HIPIO;
0414   // don't store userVariable in main, to save time
0415   if (!isCollector)
0416     HIPIO.writeHIPUserVariables(theAlignables, suvarfile.c_str(), theIteration, false, ioerr);
0417 
0418   // now calculate alignment corrections...
0419   int ialigned = 0;
0420   // iterate over alignment parameters
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       // set these parameters 'valid'
0446       ali->alignmentParameters()->setValid(true);
0447       // increase counter
0448       ialigned++;
0449     } else
0450       par->setValid(false);
0451   }
0452   //end looping over alignables
0453 
0454   edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::terminate] Aligned units: " << ialigned;
0455 
0456   // fill alignable wise root tree
0457   fillAlignablesMonitor(iSetup);
0458 
0459   edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Writing aligned parameters to file: " << theAlignables.size()
0460                                << ", for Iteration " << theIteration;
0461 
0462   // write user variables
0463   if (isCollector)
0464     HIPIO.writeHIPUserVariables(theAlignables, suvarfile.c_str(), theIteration, false, ioerr);
0465 
0466   // write new absolute positions to disk
0467   theIO.writeAlignableAbsolutePositions(theAlignables, salignedfile.c_str(), theIteration, false, ioerr);
0468 
0469   // write alignment parameters to disk
0470   //theIO.writeAlignmentParameters(theAlignables,
0471   //                 sparameterfile.c_str(),theIteration,false,ioerr);
0472 
0473   // write iteration number to file
0474   writeIterationFile(siterationfile, theIteration);
0475 
0476   // write out trees and close root file
0477   // Survey tree
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   // Alignable-wise tree is only filled once at iteration 1
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   // Eventwise and hitwise monitoring trees
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   // get trajectory impact point
0523   LocalPoint alvec = tsos.localPosition();
0524   AlgebraicVector pos(hitDim);
0525   pos[0] = alvec.x();
0526 
0527   // get impact point covariance
0528   AlgebraicSymMatrix ipcovmat(hitDim);
0529   ipcovmat[0][0] = tsos.localError().positionError().xx();
0530 
0531   // get hit local position and covariance
0532   AlgebraicVector coor(hitDim);
0533   coor[0] = hit->localPosition().x();
0534 
0535   AlgebraicSymMatrix covmat(hitDim);
0536   covmat[0][0] = hit->localPositionError().xx();
0537 
0538   // add hit and impact point covariance matrices
0539   covmat = covmat + ipcovmat;
0540 
0541   // calculate the x pull of this hit
0542   double xpull = 0.;
0543   if (covmat[0][0] != 0.)
0544     xpull = (pos[0] - coor[0]) / sqrt(fabs(covmat[0][0]));
0545 
0546   // get Alignment Parameters
0547   AlignmentParameters* params = ali->alignmentParameters();
0548   HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
0549   uservar->datatype = theDataGroup;
0550   // get derivatives
0551   AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
0552   // calculate user parameters
0553   int npar = derivs2D.num_row();
0554   AlgebraicMatrix derivs(npar, hitDim, 0);  // This is jT
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   // invert covariance matrix
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   //for alignable chi squared
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   // get trajectory impact point
0625   LocalPoint alvec = tsos.localPosition();
0626   AlgebraicVector pos(hitDim);
0627   pos[0] = alvec.x();
0628   pos[1] = alvec.y();
0629 
0630   // get impact point covariance
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   // get hit local position and covariance
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   // add hit and impact point covariance matrices
0647   covmat = covmat + ipcovmat;
0648 
0649   // calculate the x pull and y pull of this hit
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   // get Alignment Parameters
0658   AlignmentParameters* params = ali->alignmentParameters();
0659   HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
0660   uservar->datatype = theDataGroup;
0661   // get derivatives
0662   AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
0663   // calculate user parameters
0664   int npar = derivs2D.num_row();
0665   AlgebraicMatrix derivs(npar, hitDim, 0);  // This is jT
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   // invert covariance matrix
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   //for alignable chi squared
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 // Run the algorithm on trajectories and tracks -------------------------------
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   // loop over tracks
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       // Reweight by the specified eta distribution
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     // fill track parameters in root tree
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     // loop over measurements
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       // const TransientTrackingRecHit* ttrhit = &(*meas.recHit());
0807       // const TrackingRecHit *hit = ttrhit->hit();
0808       const TransientTrackingRecHit* hit = &(*meas.recHit());
0809 
0810       if (hit->isValid() && theAlignableDetAccessor->detAndSubdetInMap(hit->geographicalId())) {
0811         // this is the updated state (including the current hit)
0812         // TrajectoryStateOnSurface tsos=meas.updatedState();
0813         // combine fwd and bwd predicted state to get state
0814         // which excludes current hit
0815 
0816         //////////Hit prescaling part
0817         if (eventInfo.clusterValueMap()) {
0818           // check from the PrescalingMap if the hit was taken.
0819           // If not skip to the next TM
0820           // bool hitTaken=false;
0821           AlignmentClusterFlag myflag;
0822 
0823           int subDet = hit->geographicalId().subdetId();
0824           //take the actual RecHit out of the Transient one
0825           const TrackingRecHit* rechit = hit->hit();
0826           if (subDet > 2) {  // AM: if possible use enum instead of hard-coded value
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                 // myflag=PrescMap[stripclust];
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             }  //end if type = SiStripRecHit1D
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                 // myflag=PrescMap[stripclust];
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             }  //end if type == SiStripRecHit2D
0851           }    //end if hit from strips
0852           else {
0853             const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(rechit);
0854             if (pixelhit) {
0855               SiPixelClusterRefNew pixelclust(pixelhit->cluster());
0856               // myflag=PrescMap[pixelclust];
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           }  //end 'else' it is a pixel hit
0863           if (!myflag.isTaken())
0864             continue;
0865         }  //end if Prescaled Hits
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       }  //hit valid
0876     }
0877 
0878     // transform RecHit vector to AlignableDet vector
0879     std::vector<AlignableDetOrUnitPtr> alidetvec = theAlignableDetAccessor->alignablesFromHits(hitvec);
0880 
0881     // get concatenated alignment parameters for list of alignables
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     // loop over vectors(hit,tsos)
0888     while (itsos != tsosvec.end()) {
0889       // get AlignableDet for this hit
0890       const GeomDet* det = (*ihit)->det();
0891       // int subDet= (*ihit)->geographicalId().subdetId();
0892       uint32_t nhitDim = (*ihit)->dimension();
0893 
0894       AlignableDetOrUnitPtr alidet = theAlignableDetAccessor->alignableFromGeomDet(det);
0895 
0896       // get relevant Alignable
0897       Alignable* ali = aap.alignableFromAlignableDet(alidet);
0898 
0899       if (ali != nullptr) {
0900         const HIPAlignableSpecificParameters* alispecifics = findAlignableSpecs(ali);
0901         const TrajectoryStateOnSurface& tsos = *itsos;
0902 
0903         //  LocalVector v = tsos.localDirection();
0904         //  double proj_z = v.dot(LocalVector(0,0,1));
0905 
0906         //In fact, sin_theta=Abs(mom_z)
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         //Make cut on hit impact angle, reduce collision hits perpendicular to modules
0915         if (IsCollision) {
0916           if (angle > col_cut)
0917             alihitwt = 0;
0918         } else {
0919           if (angle < cos_cut)
0920             alihitwt = 0;
0921         }
0922 
0923         // Fill hit monitor variables
0924         theMonitorConfig.hitmonitorvars.m_angle = angle;
0925         theMonitorConfig.hitmonitorvars.m_sinTheta = sin_theta;
0926         theMonitorConfig.hitmonitorvars.m_detId = ali->id();
0927 
0928         // Check pixel XY and Q probabilities
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               // Prob X, Y are deprecated
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   }  // end of track loop
0978 
0979   // fill eventwise root tree (with prescale defined in pset)
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 // set alignment position error
1018 void HIPAlignmentAlgorithm::setAlignmentPositionError(void) {
1019   // Check if user wants to override APE
1020   if (!theApplyAPE) {
1021     edm::LogInfo("Alignment") << "[HIPAlignmentAlgorithm::setAlignmentPositionError] No APE applied";
1022     return;  // NO APE APPLIED
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     // Printout for debug
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     // set APE
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 // calculate APE
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);  // should have been caught in the constructor
1082 }
1083 
1084 // ----------------------------------------------------------------------------
1085 // book root trees
1086 void HIPAlignmentAlgorithm::bookRoot(void) {
1087   // create ROOT files
1088   if (doTrackHitMonitoring && !isCollector) {
1089     theTrackHitMonitorIORootFile = TFile::Open(theMonitorConfig.outfile.c_str(), "update");
1090     theTrackHitMonitorIORootFile->cd();
1091     // book event-wise ROOT Tree
1092     if (theMonitorConfig.fillTrackMonitoring) {
1093       TString tname = Form("T1_%i", theIteration);
1094       theTrackMonitorTree = new TTree(tname, "Eventwise tree");
1095       //theTrackMonitorTree->Branch("Run",     &m_Run,     "Run/I");
1096       //theTrackMonitorTree->Branch("Event",   &m_Event,   "Event/I");
1097       theTrackMonitorTree->Branch("DataType", &m_datatype);
1098       theMonitorConfig.trackmonitorvars.setTree(theTrackMonitorTree);
1099       theMonitorConfig.trackmonitorvars.bookBranches();
1100     }
1101     // book hit-wise ROOT Tree
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   // book alignable-wise ROOT Tree
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   // book survey-wise ROOT Tree only if survey is enabled
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 // fill alignable-wise root tree
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   //Retrieve tracker topology from geometry
1156   const TrackerTopology* const tTopo = &iSetup.getData(topoToken2_);
1157 
1158   for (const auto& ali : theAlignables) {
1159     AlignmentParameters* dap = ali->alignmentParameters();
1160 
1161     // consider only those parameters classified as 'valid'
1162     if (dap->isValid()) {
1163       // get number of hits from user variable
1164       HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(dap->userVariables());
1165       m2_Nhit = uservar->nhit;
1166       m2_datatype = uservar->datatype;
1167 
1168       // get type/layer
1169       std::pair<int, int> tl = theAlignmentParameterStore->typeAndLayer(ali, tTopo);
1170       m2_Type = tl.first;
1171       m2_Layer = tl.second;
1172 
1173       // get identifier (as for IO)
1174       m2_Id = ali->id();
1175       m2_ObjId = ali->alignableObjectId();
1176 
1177       // get position
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   // Alignment parameters
1228   AlignmentParameters* par = ali->alignmentParameters();
1229   const HIPAlignableSpecificParameters* alispecifics = findAlignableSpecs(ali);
1230   // access user variables
1231   HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(par->userVariables());
1232   int nhit = uservar->nhit;
1233   // The following variable is needed for the extended 1D/2D hit fix using
1234   // matrix shrinkage and expansion
1235   // int hitdim = uservar->hitdim;
1236 
1237   // Test nhits
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   // these are the alignment corrections+covariance (for selected params)
1253   int npar = jtve.num_row();
1254   AlgebraicVector params(npar);
1255   AlgebraicVector paramerr(npar);
1256   AlgebraicSymMatrix cov(npar);
1257 
1258   // errors of parameters
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     // Counting step for per-alignable reweighting
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();  // This vector should have 1-to-1 correspondence with the alignables vector
1356       for (const auto& ali : theAlignables) {
1357         // No need for the user variables already attached to the alignables
1358         // Just count from what you read.
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()) {  // This is a new alignable
1364             std::unordered_map<int, pawt_t> newmap;
1365             ali_datatypecountpair_map[ali] = newmap;
1366             ali_datatypecountpair_map[ali][alijobdtype] = alijobnhits;
1367           } else {  // Alignable already exists in the map
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;  // Delete new user variables as they are added
1375         }
1376         iuvar++;
1377       }  // End loop over alignables
1378     }    // End loop over subjobs
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     // add
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;  // Set the data type of alignable to that specified for the collector job (-2 by default)
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;  // Delete new user variables as they are added
1454       }
1455 
1456       uvarvecadd.push_back(uvar);
1457       iuvarnew++;
1458     }
1459 
1460     theAlignmentParameterStore->attachUserVariables(theAlignables, uvarvecadd, ioerr);
1461 
1462     // fill Eventwise Tree
1463     if (doTrackHitMonitoring) {
1464       uvfile = theCollectorPath + "/job" + str + "/" + theMonitorConfig.outfilecore;
1465       monitorFileList.push_back(uvfile);
1466     }
1467   }  // end loop on jobs
1468 
1469   // Collect monitor (eventwise and hitwise) trees
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) {  // This should never happen
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   // Leave it to HIPAlignmentAlgorithm::terminate to write the trees and close theTrackHitMonitorIORootFile
1524 
1525   delete hittrees;
1526   delete eventtrees;
1527   for (TFile*& finput : finputlist)
1528     finput->Close();
1529 
1530   // Rename the trees to standard names
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!