CommonAlignmentAlgorithm/plugins/ApeSettingAlgorithm.cc

0001 /**
0002  * \file MillePedeAlignmentAlgorithm.cc
0003  *
0004  *  \author    : Gero Flucke/Ivan Reid
0005  *  date       : February 2009 *  $Revision: 1.11 $
0006  *  $Date: 2011/03/22 09:49:50 $
0007  *  (last update by $Author: innocent $)
0008  */
0009 /*
0010  * The APE record and the ASCII file contain the covariance matrix elements
0011  * in units of cm^2
0012  *# Parameters:
0013  *#    saveApeToASCII -- Do we write out an APE text file?
0014  *#    saveComposites -- Do we write APEs for composite detectors?
0015  *#    saveLocalNotGlobal -- Do we write the APEs in the local or global coordinates?
0016  *#    apeASCIISaveFile -- The name of the save-file.
0017  *#    readApeFromASCII -- Do we read in APEs from a text file?
0018  *#    readLocalNotGlobal -- Do we read APEs in the local or the global frame?
0019  *#    readFullLocalMatrix -- Do we read the full local matrix or just the diagonal elements?
0020  *#                        -- Always write full matrix
0021  *# Full matrix format: DetID dxx dxy dyy dxz dyz dzz
0022  *# Diagonal element format: DetID sqrt(dxx) sqrt(dyy) sqrt(dzz)
0023  *#    setComposites -- Do we set the APEs for composite detectors or just ignore them?
0024  *#    apeASCIIReadFile -- Input file name.
0025  *# Also note:
0026  *#    process.AlignmentProducer.saveApeToDB -- to save as an sqlite file
0027  *# and associated entries in _cfg.py
0028  */
0029
0030 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
0031 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmBase.h"
0032 #include "Alignment/CommonAlignment/interface/AlignableModifier.h"
0033
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035
0036 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0037 #include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
0038 #include "CLHEP/Matrix/SymMatrix.h"
0039
0040 #include <fstream>
0041 #include <string>
0042 #include <set>
0043
0044 #include "DataFormats/GeometryCommonDetAlgo/interface/AlignmentPositionError.h"
0045 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
0046
0047 #include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
0048 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
0049 #include "Alignment/CommonAlignment/interface/Alignable.h"
0050
0051 // includes to make known that they inherit from Alignable:
0052 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0053 #include "Alignment/MuonAlignment/interface/AlignableMuon.h"
0054 #include "Alignment/CommonAlignment/interface/AlignableExtras.h"
0055
0056 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0057
0058 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0059 #include "FWCore/Framework/interface/ConsumesCollector.h"
0060
0061 class ApeSettingAlgorithm : public AlignmentAlgorithmBase {
0062 public:
0063   /// Constructor
0064   ApeSettingAlgorithm(const edm::ParameterSet &cfg, const edm::ConsumesCollector &iC);
0065
0066   /// Destructor
0067   ~ApeSettingAlgorithm() override;
0068
0069   /// Call at beginning of job
0070   void initialize(const edm::EventSetup &setup,
0071                   AlignableTracker *tracker,
0072                   AlignableMuon *muon,
0073                   AlignableExtras *extras,
0074                   AlignmentParameterStore *store) override;
0075
0076   /// Call at end of job
0077   void terminate(const edm::EventSetup &iSetup) override;
0078
0079   /// Run the algorithm
0080   void run(const edm::EventSetup &setup, const EventInfo &eventInfo) override;
0081
0082 private:
0083   edm::ParameterSet theConfig;
0084   AlignableNavigator *theAlignableNavigator;
0085   AlignableTracker *theTracker;
0086   bool saveApeToAscii_, readApeFromAscii_, readFullLocalMatrix_;
0087   bool readLocalNotGlobal_, saveLocalNotGlobal_;
0088   bool setComposites_, saveComposites_;
0089 };
0090
0091 //____________________________________________________
0092 //____________________________________________________
0093 //____________________________________________________
0094 //____________________________________________________
0095
0096 // Constructor ----------------------------------------------------------------
0097 //____________________________________________________
0098 ApeSettingAlgorithm::ApeSettingAlgorithm(const edm::ParameterSet &cfg, const edm::ConsumesCollector &iC)
0099     : AlignmentAlgorithmBase(cfg, iC), theConfig(cfg), theAlignableNavigator(nullptr) {
0100   edm::LogInfo("Alignment") << "@SUB=ApeSettingAlgorithm"
0101                             << "Start.";
0102   saveApeToAscii_ = theConfig.getUntrackedParameter<bool>("saveApeToASCII");
0103   saveComposites_ = theConfig.getUntrackedParameter<bool>("saveComposites");
0104   saveLocalNotGlobal_ = theConfig.getUntrackedParameter<bool>("saveLocalNotGlobal");
0105   readApeFromAscii_ = theConfig.getParameter<bool>("readApeFromASCII");
0106   readLocalNotGlobal_ = theConfig.getParameter<bool>("readLocalNotGlobal");
0107   readFullLocalMatrix_ = theConfig.getParameter<bool>("readFullLocalMatrix");
0108   setComposites_ = theConfig.getParameter<bool>("setComposites");
0109 }
0110
0111 // Destructor ----------------------------------------------------------------
0112 //____________________________________________________
0113 ApeSettingAlgorithm::~ApeSettingAlgorithm() { delete theAlignableNavigator; }
0114
0115 // Call at beginning of job ---------------------------------------------------
0116 //____________________________________________________
0117 void ApeSettingAlgorithm::initialize(const edm::EventSetup &setup,
0118                                      AlignableTracker *tracker,
0119                                      AlignableMuon *muon,
0120                                      AlignableExtras *extras,
0121                                      AlignmentParameterStore *store) {
0122   theAlignableNavigator = new AlignableNavigator(tracker, muon);
0123   theTracker = tracker;
0124
0125   if (readApeFromAscii_) {
0126     std::ifstream apeReadFile(
0127         theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath().c_str());  //requires <fstream>
0128     if (!apeReadFile.good()) {
0129       edm::LogInfo("Alignment") << "@SUB=initialize"
0130                                 << "Problem opening APE file: skipping"
0131                                 << theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath();
0132       return;
0133     }
0134     std::set<int> apeList;  //To avoid duplicates
0135     while (!apeReadFile.eof()) {
0136       int apeId = 0;
0137       double x11, x21, x22, x31, x32, x33, ignore;
0138       if (!readLocalNotGlobal_ || readFullLocalMatrix_) {
0139         apeReadFile >> apeId >> x11 >> x21 >> x22 >> x31 >> x32 >> x33 >> ignore >> ignore >> ignore >> ignore >>
0140             ignore >> ignore >> ignore >> ignore >> ignore >> ignore >> ignore >> ignore >> ignore >> ignore >>
0141             ignore >> std::ws;
0142       } else {
0143         apeReadFile >> apeId >> x11 >> x22 >> x33 >> std::ws;
0144       }
0145       //idr What sanity checks do we need to put here?
0146       if (apeId != 0)  //read appears valid?
0147       {
0148         if (apeList.find(apeId) == apeList.end())  //Not previously done
0149         {
0150           DetId id(apeId);
0151           AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(id));  //NULL if none
0152           if (alidet) {
0153             if ((alidet->components().empty()) || setComposites_)  //the problem with glued dets...
0154             {
0155               GlobalErrorExtended globErr;
0156               if (readLocalNotGlobal_) {
0157                 AlgebraicSymMatrix33 as;
0158                 if (readFullLocalMatrix_) {
0159                   as[0][0] = x11;
0160                   as[1][0] = x21;
0161                   as[1][1] = x22;
0162                   as[2][0] = x31;
0163                   as[2][1] = x32;
0164                   as[2][2] = x33;
0165                 } else {
0166                   as[0][0] = x11 * x11;
0167                   as[1][1] = x22 * x22;
0168                   as[2][2] = x33 * x33;
0169                 }  //local cov.
0170                 align::RotationType rt = alidet->globalRotation();
0171                 AlgebraicMatrix33 am;
0172                 am[0][0] = rt.xx();
0173                 am[0][1] = rt.xy();
0174                 am[0][2] = rt.xz();
0175                 am[1][0] = rt.yx();
0176                 am[1][1] = rt.yy();
0177                 am[1][2] = rt.yz();
0178                 am[2][0] = rt.zx();
0179                 am[2][1] = rt.zy();
0180                 am[2][2] = rt.zz();
0181                 globErr = GlobalErrorExtended(ROOT::Math::SimilarityT(am, as));
0182               } else {
0183                 if (readFullLocalMatrix_)
0184                   globErr =
0185                       GlobalErrorExtended(x11, x21, x31, 0, 0, 0, x22, x32, 0, 0, 0, x33, 0, 0, 0, 0, 0, 0, 0, 0, 0);
0186                 else {
0187                   globErr = GlobalErrorExtended(
0188                       x11 * x11, 0, 0, 0, 0, 0, x22 * x22, 0, 0, 0, 0, x33 * x33, 0, 0, 0, 0, 0, 0, 0, 0, 0);
0189                 }
0190               }
0191               alidet->setAlignmentPositionError(globErr, false);  // do not propagate down!
0192               apeList.insert(apeId);                              //Flag it's been set
0193             } else {
0194               edm::LogInfo("Alignment") << "@SUB=initialize"
0195                                         << "Not Setting APE for Composite DetId " << apeId;
0196             }
0197           }
0198         } else {
0199           edm::LogInfo("Alignment") << "@SUB=initialize"
0200                                     << "Skipping duplicate APE for DetId " << apeId;
0201         }
0202       }
0203     }
0204     apeReadFile.close();
0205     edm::LogInfo("Alignment") << "@SUB=initialize"
0206                               << "Set " << apeList.size() << " APE values.";
0207   }
0208 }
0209
0210 // Call at end of job ---------------------------------------------------------
0211 //____________________________________________________
0212 void ApeSettingAlgorithm::terminate(const edm::EventSetup &iSetup) {
0213   if (saveApeToAscii_) {
0214     AlignmentErrorsExtended *aliErr = theTracker->alignmentErrors();
0215     int theSize = aliErr->m_alignError.size();
0216     std::ofstream apeSaveFile(
0217         theConfig.getUntrackedParameter<std::string>("apeASCIISaveFile").c_str());  //requires <fstream>
0218     for (int i = 0; i < theSize; ++i) {
0219       int id = aliErr->m_alignError[i].rawId();
0220       AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(DetId(id)));  //NULL if none
0221       if (alidet && ((alidet->components().empty()) || saveComposites_)) {
0222         apeSaveFile << id;
0223         CLHEP::HepSymMatrix sm = aliErr->m_alignError[i].matrix();
0224         if (saveLocalNotGlobal_) {
0225           align::RotationType rt = alidet->globalRotation();
0226           AlgebraicMatrix am(3, 3);
0227           am[0][0] = rt.xx();
0228           am[0][1] = rt.xy();
0229           am[0][2] = rt.xz();
0230           am[1][0] = rt.yx();
0231           am[1][1] = rt.yy();
0232           am[1][2] = rt.yz();
0233           am[2][0] = rt.zx();
0234           am[2][1] = rt.zy();
0235           am[2][2] = rt.zz();
0236           sm = sm.similarity(am);  //symmetric matrix
0237         }  //transform to local
0238         for (int j = 0; j < sm.num_row(); ++j)
0239           for (int k = 0; k <= j; ++k)
0240             apeSaveFile << "  " << sm[j][k];  //always write full matrix
0241
0242         apeSaveFile << std::endl;
0243       }
0244     }
0245     delete aliErr;
0246     apeSaveFile.close();
0247   }
0248   // clean up at end:  // FIXME: should we delete here or in destructor?
0249   delete theAlignableNavigator;
0250   theAlignableNavigator = nullptr;
0251 }
0252
0253 // Run the algorithm on trajectories and tracks -------------------------------
0254 //____________________________________________________
0255 void ApeSettingAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo) {
0256   // nothing to do here?
0257 }
0258
0259 // Plugin definition for the algorithm
0260 DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory, ApeSettingAlgorithm, "ApeSettingAlgorithm");