Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 11:57:09

0001 #include "CondFormats/Alignment/interface/Alignments.h"
0002 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0003 #include "CondFormats/Alignment/interface/AlignmentSurfaceDeformations.h"
0004 #include "CondFormats/Alignment/interface/Definitions.h"
0005 #include "CLHEP/Vector/RotationInterfaces.h"
0006 #include "CondFormats/AlignmentRecord/interface/TrackerSurveyRcd.h"
0007 #include "CondFormats/AlignmentRecord/interface/TrackerSurveyErrorExtendedRcd.h"
0008 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0009 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESTransientHandle.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 
0018 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0019 #include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
0020 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0021 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0022 #include "CondCore/SiPixelPlugins/interface/SiPixelPayloadInspectorHelper.h"
0023 
0024 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
0025 #include "Geometry/CommonTopologies/interface/SurfaceDeformationFactory.h"
0026 #include "Geometry/CommonTopologies/interface/SurfaceDeformation.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
0029 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0030 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0031 #include "Geometry/CommonTopologies/interface/GeometryAligner.h"
0032 #include "Alignment/CommonAlignment/interface/Utilities.h"
0033 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
0034 #include "Alignment/CommonAlignment/interface/Alignable.h"
0035 #include "DataFormats/DetId/interface/DetId.h"
0036 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0037 #include "CondFormats/Alignment/interface/DetectorGlobalPosition.h"
0038 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0039 
0040 #include "TrackerGeometryCompare.h"
0041 #include "TFile.h"
0042 #include "CLHEP/Vector/ThreeVector.h"
0043 
0044 // Database
0045 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0046 #include "FWCore/ServiceRegistry/interface/Service.h"
0047 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0048 
0049 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0050 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0051 
0052 #include <iostream>
0053 #include <fstream>
0054 #include <sstream>
0055 
0056 TrackerGeometryCompare::TrackerGeometryCompare(const edm::ParameterSet& cfg)
0057     : cpvTokenDDD_(esConsumes()),
0058       cpvTokenDD4hep_(esConsumes()),
0059       topoToken_(esConsumes()),
0060       geomDetToken_(esConsumes()),
0061       ptpToken_(esConsumes()),
0062       pixQualityToken_(esConsumes()),
0063       stripQualityToken_(esConsumes()),
0064       referenceTracker(nullptr),
0065       dummyTracker(nullptr),
0066       currentTracker(nullptr),
0067       theSurveyIndex(0),
0068       theSurveyValues(nullptr),
0069       theSurveyErrors(nullptr),
0070       levelStrings_(cfg.getUntrackedParameter<std::vector<std::string> >("levels")),
0071       fromDD4hep_(cfg.getUntrackedParameter<bool>("fromDD4hep")),
0072       writeToDB_(cfg.getUntrackedParameter<bool>("writeToDB")),
0073       commonTrackerLevel_(align::invalid),
0074       moduleListFile_(),
0075       moduleList_(0),
0076       inputRootFile1_(nullptr),
0077       inputRootFile2_(nullptr),
0078       inputTree01_(nullptr),
0079       inputTree02_(nullptr),
0080       inputTree11_(nullptr),
0081       inputTree12_(nullptr),
0082       m_nBins_(10000),
0083       m_rangeLow_(-.1),
0084       m_rangeHigh_(.1),
0085       firstEvent_(true),
0086       m_vtkmap_(13) {
0087   moduleListName_ = cfg.getUntrackedParameter<std::string>("moduleList");
0088 
0089   //input is ROOT
0090   inputFilename1_ = cfg.getUntrackedParameter<std::string>("inputROOTFile1");
0091   inputFilename2_ = cfg.getUntrackedParameter<std::string>("inputROOTFile2");
0092   inputTreenameAlign_ = cfg.getUntrackedParameter<std::string>("treeNameAlign");
0093   inputTreenameDeform_ = cfg.getUntrackedParameter<std::string>("treeNameDeform");
0094 
0095   //output file
0096   filename_ = cfg.getUntrackedParameter<std::string>("outputFile");
0097   //output dir for surface deformations
0098   surfdir_ = cfg.getUntrackedParameter<std::string>("surfDir");
0099 
0100   weightBy_ = cfg.getUntrackedParameter<std::string>("weightBy");
0101   setCommonTrackerSystem_ = cfg.getUntrackedParameter<std::string>("setCommonTrackerSystem");
0102   detIdFlag_ = cfg.getUntrackedParameter<bool>("detIdFlag");
0103   detIdFlagFile_ = cfg.getUntrackedParameter<std::string>("detIdFlagFile");
0104   weightById_ = cfg.getUntrackedParameter<bool>("weightById");
0105   weightByIdFile_ = cfg.getUntrackedParameter<std::string>("weightByIdFile");
0106 
0107   // if want to use, make id cut list
0108   if (detIdFlag_) {
0109     std::ifstream fin;
0110     fin.open(detIdFlagFile_.c_str());
0111 
0112     while (!fin.eof() && fin.good()) {
0113       uint32_t id;
0114       fin >> id;
0115       detIdFlagVector_.push_back(id);
0116     }
0117     fin.close();
0118   }
0119 
0120   // turn weightByIdFile into weightByIdVector
0121   if (weightById_) {
0122     std::ifstream inFile;
0123     inFile.open(weightByIdFile_.c_str());
0124     while (!inFile.eof()) {
0125       unsigned int listId;
0126       inFile >> listId;
0127       inFile.ignore(256, '\n');
0128 
0129       weightByIdVector_.push_back(listId);
0130     }
0131     inFile.close();
0132   }
0133 
0134   //root configuration
0135   theFile_ = new TFile(filename_.c_str(), "RECREATE");
0136   alignTree_ = new TTree("alignTree",
0137                          "alignTree");  //,"id:level:mid:mlevel:sublevel:x:y:z:r:phi:a:b:c:dx:dy:dz:dr:dphi:da:db:dc");
0138   alignTree_->Branch("id", &id_, "id/I");
0139   alignTree_->Branch("badModuleQuality", &badModuleQuality_, "badModuleQuality/I");
0140   alignTree_->Branch("inModuleList", &inModuleList_, "inModuleList/I");
0141   alignTree_->Branch("level", &level_, "level/I");
0142   alignTree_->Branch("mid", &mid_, "mid/I");
0143   alignTree_->Branch("mlevel", &mlevel_, "mlevel/I");
0144   alignTree_->Branch("sublevel", &sublevel_, "sublevel/I");
0145   alignTree_->Branch("x", &xVal_, "x/F");
0146   alignTree_->Branch("y", &yVal_, "y/F");
0147   alignTree_->Branch("z", &zVal_, "z/F");
0148   alignTree_->Branch("r", &rVal_, "r/F");
0149   alignTree_->Branch("phi", &phiVal_, "phi/F");
0150   alignTree_->Branch("eta", &etaVal_, "eta/F");
0151   alignTree_->Branch("alpha", &alphaVal_, "alpha/F");
0152   alignTree_->Branch("beta", &betaVal_, "beta/F");
0153   alignTree_->Branch("gamma", &gammaVal_, "gamma/F");
0154   alignTree_->Branch("dx", &dxVal_, "dx/F");
0155   alignTree_->Branch("dy", &dyVal_, "dy/F");
0156   alignTree_->Branch("dz", &dzVal_, "dz/F");
0157   alignTree_->Branch("dr", &drVal_, "dr/F");
0158   alignTree_->Branch("dphi", &dphiVal_, "dphi/F");
0159   alignTree_->Branch("dalpha", &dalphaVal_, "dalpha/F");
0160   alignTree_->Branch("dbeta", &dbetaVal_, "dbeta/F");
0161   alignTree_->Branch("dgamma", &dgammaVal_, "dgamma/F");
0162   alignTree_->Branch("du", &duVal_, "du/F");
0163   alignTree_->Branch("dv", &dvVal_, "dv/F");
0164   alignTree_->Branch("dw", &dwVal_, "dw/F");
0165   alignTree_->Branch("da", &daVal_, "da/F");
0166   alignTree_->Branch("db", &dbVal_, "db/F");
0167   alignTree_->Branch("dg", &dgVal_, "dg/F");
0168   alignTree_->Branch("useDetId", &useDetId_, "useDetId/I");
0169   alignTree_->Branch("detDim", &detDim_, "detDim/I");
0170   alignTree_->Branch("surW", &surWidth_, "surW/F");
0171   alignTree_->Branch("surL", &surLength_, "surL/F");
0172   alignTree_->Branch("surRot", &surRot_, "surRot[9]/D");
0173   alignTree_->Branch("identifiers", &identifiers_, "identifiers[6]/I");
0174   alignTree_->Branch("type", &type_, "type/I");
0175   alignTree_->Branch("surfDeform", &surfDeform_, "surfDeform[13]/D");
0176 
0177   for (std::vector<TrackerMap>::iterator it = m_vtkmap_.begin(); it != m_vtkmap_.end(); ++it) {
0178     it->setPalette(1);
0179     it->addPixel(true);
0180   }
0181 
0182   edm::Service<TFileService> fs;
0183   TFileDirectory subDir_All = fs->mkdir("AllSubdetectors");
0184   TFileDirectory subDir_PXB = fs->mkdir("PixelBarrel");
0185   TFileDirectory subDir_PXF = fs->mkdir("PixelEndcap");
0186   for (int ii = 0; ii < 13; ++ii) {
0187     std::stringstream histname0;
0188     histname0 << "SurfDeform_Par_" << ii;
0189     m_h1_[histname0.str()] = subDir_All.make<TH1D>(
0190         (histname0.str()).c_str(), (histname0.str()).c_str(), m_nBins_, m_rangeLow_, m_rangeHigh_);
0191 
0192     std::stringstream histname1;
0193     histname1 << "SurfDeform_PixelBarrel_Par_" << ii;
0194     m_h1_[histname1.str()] = subDir_PXB.make<TH1D>(
0195         (histname1.str()).c_str(), (histname1.str()).c_str(), m_nBins_, m_rangeLow_, m_rangeHigh_);
0196 
0197     std::stringstream histname2;
0198     histname2 << "SurfDeform_PixelEndcap_Par_" << ii;
0199     m_h1_[histname2.str()] = subDir_PXF.make<TH1D>(
0200         (histname2.str()).c_str(), (histname2.str()).c_str(), m_nBins_, m_rangeLow_, m_rangeHigh_);
0201   }
0202 }
0203 
0204 void TrackerGeometryCompare::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0205   edm::ParameterSetDescription desc;
0206   desc.setComment("Validates alignment payloads by comparing positions of tracker modules positiona and orientations");
0207   desc.addUntracked<std::vector<std::string> >("levels", {});
0208   desc.addUntracked<bool>("fromDD4hep", false);
0209   desc.addUntracked<bool>("writeToDB", false);
0210   desc.addUntracked<std::string>("moduleList", "moduleList.txt");
0211   desc.addUntracked<std::string>("inputROOTFile1", "IDEAL");
0212   desc.addUntracked<std::string>("inputROOTFile2", "idealtracker2.root");
0213   desc.addUntracked<std::string>("treeNameAlign", "alignTree");
0214   desc.addUntracked<std::string>("treeNameDeform", "alignTreeDeformations");
0215   desc.addUntracked<std::string>("outputFile", "output.root");
0216   desc.addUntracked<std::string>("surfDir", ".");
0217   desc.addUntracked<std::string>("weightBy", "DetUnit");
0218   desc.addUntracked<std::string>("setCommonTrackerSystem", "NONE");
0219   desc.addUntracked<bool>("detIdFlag", false);
0220   desc.addUntracked<std::string>("detIdFlagFile", "blah.txt");
0221   desc.addUntracked<bool>("weightById", false);
0222   desc.addUntracked<std::string>("weightByIdFile", "blah2.txt");
0223   descriptions.addWithDefaultLabel(desc);
0224 }
0225 
0226 void TrackerGeometryCompare::beginJob() { firstEvent_ = true; }
0227 
0228 void TrackerGeometryCompare::endJob() {
0229   int iname(0);
0230   for (std::vector<TrackerMap>::iterator it = m_vtkmap_.begin(); it != m_vtkmap_.end(); ++it) {
0231     std::stringstream mapname;
0232     mapname << surfdir_ << "/TkMap_SurfDeform_0" << iname << ".png";
0233     it->save(true, 0, 0, mapname.str());
0234     mapname.str(std::string());
0235     mapname.clear();
0236     mapname << surfdir_ << "/TkMap_SurfDeform_1" << iname << ".pdf";
0237     it->save(true, 0, 0, mapname.str());
0238     ++iname;
0239   }
0240 
0241   theFile_->cd();
0242   alignTree_->Write();
0243   theFile_->Close();
0244 }
0245 
0246 void TrackerGeometryCompare::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0247   if (firstEvent_) {
0248     //Retrieve tracker topology from geometry
0249     const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0250 
0251     //upload the ROOT geometries
0252     createROOTGeometry(iSetup);
0253 
0254     //setting the levels being used in the geometry comparator
0255     edm::LogInfo("TrackerGeometryCompare") << "levels: " << levelStrings_.size();
0256     for (const auto& level : levelStrings_) {
0257       m_theLevels.push_back(currentTracker->objectIdProvider().stringToId(level));
0258       edm::LogInfo("TrackerGeometryCompare") << "level: " << level;
0259       edm::LogInfo("TrackerGeometryCompare")
0260           << "structure type: " << currentTracker->objectIdProvider().stringToId(level);
0261     }
0262 
0263     //set common tracker system first
0264     // if setting the tracker common system
0265     if (setCommonTrackerSystem_ != "NONE") {
0266       setCommonTrackerSystem();
0267     }
0268 
0269     //compare the goemetries
0270     compareGeometries(referenceTracker, currentTracker, tTopo, iSetup);
0271     compareSurfaceDeformations(inputTree11_, inputTree12_);
0272 
0273     //write out ntuple
0274     //might be better to do within output module
0275 
0276     if (writeToDB_) {
0277       Alignments* myAlignments = currentTracker->alignments();
0278       AlignmentErrorsExtended* myAlignmentErrorsExtended = currentTracker->alignmentErrors();
0279 
0280       // 2. Store alignment[Error]s to DB
0281       edm::Service<cond::service::PoolDBOutputService> poolDbService;
0282       // Call service
0283       if (!poolDbService.isAvailable())  // Die if not available
0284         throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
0285 
0286       poolDbService->writeOneIOV<Alignments>(*myAlignments, poolDbService->beginOfTime(), "TrackerAlignmentRcd");
0287       poolDbService->writeOneIOV<AlignmentErrorsExtended>(
0288           *myAlignmentErrorsExtended, poolDbService->beginOfTime(), "TrackerAlignmentErrorExtendedRcd");
0289     }
0290 
0291     firstEvent_ = false;
0292   }
0293 }
0294 
0295 void TrackerGeometryCompare::createROOTGeometry(const edm::EventSetup& iSetup) {
0296   int inputRawId1, inputRawId2;
0297   double inputX1, inputY1, inputZ1, inputX2, inputY2, inputZ2;
0298   double inputAlpha1, inputBeta1, inputGamma1, inputAlpha2, inputBeta2, inputGamma2;
0299 
0300   //Retrieve tracker topology from geometry
0301   const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0302 
0303   // Fill module IDs from file into a list
0304   moduleListFile_.open(moduleListName_);
0305   if (moduleListFile_.is_open()) {
0306     std::string line;
0307     while (!moduleListFile_.eof()) {
0308       std::getline(moduleListFile_, line);
0309       moduleList_.push_back(std::atoi(line.c_str()));
0310     }
0311   } else {
0312     edm::LogInfo("TrackerGeometryCompare") << "Error: Module list not found! Please verify that given list exists!";
0313   }
0314 
0315   //declare alignments
0316   Alignments* alignments1 = new Alignments();
0317   AlignmentErrorsExtended* alignmentErrors1 = new AlignmentErrorsExtended();
0318   if (inputFilename1_ != "IDEAL") {
0319     inputRootFile1_ = new TFile(inputFilename1_.c_str());
0320     TTree* inputTree01_ = (TTree*)inputRootFile1_->Get(inputTreenameAlign_.c_str());
0321     inputTree01_->SetBranchAddress("rawid", &inputRawId1);
0322     inputTree01_->SetBranchAddress("x", &inputX1);
0323     inputTree01_->SetBranchAddress("y", &inputY1);
0324     inputTree01_->SetBranchAddress("z", &inputZ1);
0325     inputTree01_->SetBranchAddress("alpha", &inputAlpha1);
0326     inputTree01_->SetBranchAddress("beta", &inputBeta1);
0327     inputTree01_->SetBranchAddress("gamma", &inputGamma1);
0328 
0329     int nEntries1 = inputTree01_->GetEntries();
0330     //fill alignments
0331     for (int i = 0; i < nEntries1; ++i) {
0332       inputTree01_->GetEntry(i);
0333       CLHEP::Hep3Vector translation1(inputX1, inputY1, inputZ1);
0334       CLHEP::HepEulerAngles eulerangles1(inputAlpha1, inputBeta1, inputGamma1);
0335       uint32_t detid1 = inputRawId1;
0336       AlignTransform transform1(translation1, eulerangles1, detid1);
0337       alignments1->m_align.push_back(transform1);
0338 
0339       // dummy errors
0340       // APE matrix was 3x3, now it's 6x6 (because of muons), see PR #6483
0341       CLHEP::HepSymMatrix clhepSymMatrix(6, 0);
0342       AlignTransformErrorExtended transformError(clhepSymMatrix, detid1);
0343       alignmentErrors1->m_alignError.push_back(transformError);
0344     }
0345 
0346     // to get the right order, sort by rawId
0347     std::sort(alignments1->m_align.begin(), alignments1->m_align.end());
0348     std::sort(alignmentErrors1->m_alignError.begin(), alignmentErrors1->m_alignError.end());
0349   }
0350   //------------------
0351   Alignments* alignments2 = new Alignments();
0352   AlignmentErrorsExtended* alignmentErrors2 = new AlignmentErrorsExtended();
0353   if (inputFilename2_ != "IDEAL") {
0354     inputRootFile2_ = new TFile(inputFilename2_.c_str());
0355     TTree* inputTree02_ = (TTree*)inputRootFile2_->Get(inputTreenameAlign_.c_str());
0356     inputTree02_->SetBranchAddress("rawid", &inputRawId2);
0357     inputTree02_->SetBranchAddress("x", &inputX2);
0358     inputTree02_->SetBranchAddress("y", &inputY2);
0359     inputTree02_->SetBranchAddress("z", &inputZ2);
0360     inputTree02_->SetBranchAddress("alpha", &inputAlpha2);
0361     inputTree02_->SetBranchAddress("beta", &inputBeta2);
0362     inputTree02_->SetBranchAddress("gamma", &inputGamma2);
0363 
0364     int nEntries2 = inputTree02_->GetEntries();
0365     //fill alignments
0366     for (int i = 0; i < nEntries2; ++i) {
0367       inputTree02_->GetEntry(i);
0368       CLHEP::Hep3Vector translation2(inputX2, inputY2, inputZ2);
0369       CLHEP::HepEulerAngles eulerangles2(inputAlpha2, inputBeta2, inputGamma2);
0370       uint32_t detid2 = inputRawId2;
0371       AlignTransform transform2(translation2, eulerangles2, detid2);
0372       alignments2->m_align.push_back(transform2);
0373 
0374       // dummy errors
0375       // APE matrix was 3x3, now it's 6x6 (because of muons), see PR #6483
0376       CLHEP::HepSymMatrix clhepSymMatrix(6, 0);
0377       AlignTransformErrorExtended transformError(clhepSymMatrix, detid2);
0378       alignmentErrors2->m_alignError.push_back(transformError);
0379     }
0380 
0381     // to get the right order, sort by rawId
0382     std::sort(alignments2->m_align.begin(), alignments2->m_align.end());
0383     std::sort(alignmentErrors2->m_alignError.begin(), alignmentErrors2->m_alignError.end());
0384   }
0385 
0386   //accessing the initial geometry
0387   if (!fromDD4hep_) {
0388     edm::ESTransientHandle<DDCompactView> cpv = iSetup.getTransientHandle(cpvTokenDDD_);
0389   } else {
0390     edm::ESTransientHandle<cms::DDCompactView> cpv = iSetup.getTransientHandle(cpvTokenDD4hep_);
0391   }
0392 
0393   const GeometricDet* theGeometricDet = &iSetup.getData(geomDetToken_);
0394   const PTrackerParameters* ptp = &iSetup.getData(ptpToken_);
0395   TrackerGeomBuilderFromGeometricDet trackerBuilder;
0396 
0397   //reference tracker
0398   TrackerGeometry* theRefTracker = trackerBuilder.build(theGeometricDet, *ptp, tTopo);
0399   if (inputFilename1_ != "IDEAL") {
0400     GeometryAligner aligner1;
0401     aligner1.applyAlignments<TrackerGeometry>(
0402         &(*theRefTracker), &(*alignments1), &(*alignmentErrors1), AlignTransform());
0403   }
0404   referenceTracker = new AlignableTracker(&(*theRefTracker), tTopo);
0405   //referenceTracker->setSurfaceDeformation(surfDef1, true) ;
0406 
0407   //set tracker pixel Phase, assume current tracker is the same Phase
0408   if (theRefTracker->isThere(GeomDetEnumerators::P2PXB) || theRefTracker->isThere(GeomDetEnumerators::P2PXEC)) {
0409     phase_ = SiPixelPI::phase::two;
0410   } else if (theRefTracker->isThere(GeomDetEnumerators::P1PXB) || theRefTracker->isThere(GeomDetEnumerators::P1PXEC)) {
0411     phase_ = SiPixelPI::phase::one;
0412   } else {
0413     phase_ = SiPixelPI::phase::zero;
0414   }
0415 
0416   int inputRawid1;
0417   int inputRawid2;
0418   int inputDtype1, inputDtype2;
0419   std::vector<double> inputDpar1;
0420   std::vector<double> inputDpar2;
0421   std::vector<double>* p_inputDpar1 = &inputDpar1;
0422   std::vector<double>* p_inputDpar2 = &inputDpar2;
0423 
0424   const auto& comp1 = referenceTracker->deepComponents();
0425 
0426   SurfaceDeformation* surfDef1;
0427   if (inputFilename1_ != "IDEAL") {
0428     TTree* inputTree11_ = (TTree*)inputRootFile1_->Get(inputTreenameDeform_.c_str());
0429     inputTree11_->SetBranchAddress("irawid", &inputRawid1);
0430     inputTree11_->SetBranchAddress("dtype", &inputDtype1);
0431     inputTree11_->SetBranchAddress("dpar", &p_inputDpar1);
0432 
0433     unsigned int nEntries11 = inputTree11_->GetEntries();
0434     edm::LogInfo("TrackerGeometryCompare") << " nentries11 = " << nEntries11;
0435     for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
0436       inputTree11_->GetEntry(iEntry);
0437 
0438       surfDef1 = SurfaceDeformationFactory::create(inputDtype1, inputDpar1);
0439 
0440       if (int(comp1[iEntry]->id()) == inputRawid1) {
0441         comp1[iEntry]->setSurfaceDeformation(surfDef1, true);
0442       }
0443     }
0444   }
0445 
0446   //currernt tracker
0447   TrackerGeometry* theCurTracker = trackerBuilder.build(&*theGeometricDet, *ptp, tTopo);
0448   if (inputFilename2_ != "IDEAL") {
0449     GeometryAligner aligner2;
0450     aligner2.applyAlignments<TrackerGeometry>(
0451         &(*theCurTracker), &(*alignments2), &(*alignmentErrors2), AlignTransform());
0452   }
0453   currentTracker = new AlignableTracker(&(*theCurTracker), tTopo);
0454 
0455   const auto& comp2 = currentTracker->deepComponents();
0456 
0457   SurfaceDeformation* surfDef2;
0458   if (inputFilename2_ != "IDEAL") {
0459     TTree* inputTree12_ = (TTree*)inputRootFile2_->Get(inputTreenameDeform_.c_str());
0460     inputTree12_->SetBranchAddress("irawid", &inputRawid2);
0461     inputTree12_->SetBranchAddress("dtype", &inputDtype2);
0462     inputTree12_->SetBranchAddress("dpar", &p_inputDpar2);
0463 
0464     unsigned int nEntries12 = inputTree12_->GetEntries();
0465     edm::LogInfo("TrackerGeometryCompare") << " nentries12 = " << nEntries12;
0466     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
0467       inputTree12_->GetEntry(iEntry);
0468 
0469       surfDef2 = SurfaceDeformationFactory::create(inputDtype2, inputDpar2);
0470 
0471       if (int(comp2[iEntry]->id()) == inputRawid2) {
0472         comp2[iEntry]->setSurfaceDeformation(surfDef2, true);
0473       }
0474     }
0475   }
0476 
0477   delete alignments1;
0478   delete alignmentErrors1;
0479   delete alignments2;
0480   delete alignmentErrors2;
0481 }
0482 
0483 void TrackerGeometryCompare::compareSurfaceDeformations(TTree* refTree, TTree* curTree) {
0484   if (inputFilename1_ != "IDEAL" && inputFilename2_ != "IDEAL") {
0485     int inputRawid1;
0486     int inputRawid2;
0487     int inputSubdetid1, inputSubdetid2;
0488     int inputDtype1, inputDtype2;
0489     std::vector<double> inputDpar1;
0490     std::vector<double> inputDpar2;
0491     std::vector<double>* p_inputDpar1 = &inputDpar1;
0492     std::vector<double>* p_inputDpar2 = &inputDpar2;
0493 
0494     TTree* refTree = (TTree*)inputRootFile1_->Get(inputTreenameDeform_.c_str());
0495     refTree->SetBranchAddress("irawid", &inputRawid1);
0496     refTree->SetBranchAddress("subdetid", &inputSubdetid1);
0497     refTree->SetBranchAddress("dtype", &inputDtype1);
0498     refTree->SetBranchAddress("dpar", &p_inputDpar1);
0499 
0500     TTree* curTree = (TTree*)inputRootFile2_->Get(inputTreenameDeform_.c_str());
0501     curTree->SetBranchAddress("irawid", &inputRawid2);
0502     curTree->SetBranchAddress("subdetid", &inputSubdetid2);
0503     curTree->SetBranchAddress("dtype", &inputDtype2);
0504     curTree->SetBranchAddress("dpar", &p_inputDpar2);
0505 
0506     unsigned int nEntries11 = refTree->GetEntries();
0507     unsigned int nEntries12 = curTree->GetEntries();
0508 
0509     if (nEntries11 != nEntries12) {
0510       edm::LogError("TrackerGeometryCompare") << " Surface deformation parameters in two geometries differ!\n";
0511       return;
0512     }
0513 
0514     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
0515       refTree->GetEntry(iEntry);
0516       curTree->GetEntry(iEntry);
0517       for (int ii = 0; ii < 13; ++ii) {
0518         surfDeform_[ii] = -1.0;
0519       }
0520       for (int npar = 0; npar < int(inputDpar2.size()); ++npar) {
0521         if (inputRawid1 == inputRawid2) {
0522           surfDeform_[npar] = inputDpar2.at(npar) - inputDpar1.at(npar);
0523           std::stringstream histname0;
0524           histname0 << "SurfDeform_Par_" << npar;
0525           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0526             m_h1_[histname0.str()]->Fill(surfDeform_[npar]);
0527           if (inputSubdetid1 == 1 && inputSubdetid2 == 1) {
0528             std::stringstream histname1;
0529             histname1 << "SurfDeform_PixelBarrel_Par_" << npar;
0530             if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0531               m_h1_[histname1.str()]->Fill(surfDeform_[npar]);
0532           }
0533           if (inputSubdetid1 == 2 && inputSubdetid2 == 2) {
0534             std::stringstream histname2;
0535             histname2 << "SurfDeform_PixelEndcap_Par_" << npar;
0536             if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0537               m_h1_[histname2.str()]->Fill(surfDeform_[npar]);
0538           }
0539           if (phase_ < SiPixelPI::phase::two) {
0540             if (phase_ == SiPixelPI::phase::zero) {
0541               (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0542             } else {
0543               // old-style tracker map does not support Phase1 Pixel
0544               if (!(inputSubdetid1 == 1 && inputSubdetid2 == 1) && !(inputSubdetid1 == 2 && inputSubdetid2 == 2)) {
0545                 (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0546               }
0547             }
0548           }
0549         }
0550       }
0551     }
0552 
0553   } else if (inputFilename1_ == "IDEAL" && inputFilename2_ != "IDEAL") {
0554     int inputRawid2;
0555     int inputSubdetid2;
0556     int inputDtype2;
0557     std::vector<double> inputDpar2;
0558     std::vector<double>* p_inputDpar2 = &inputDpar2;
0559 
0560     TTree* curTree = (TTree*)inputRootFile2_->Get(inputTreenameDeform_.c_str());
0561     curTree->SetBranchAddress("irawid", &inputRawid2);
0562     curTree->SetBranchAddress("subdetid", &inputSubdetid2);
0563     curTree->SetBranchAddress("dtype", &inputDtype2);
0564     curTree->SetBranchAddress("dpar", &p_inputDpar2);
0565 
0566     unsigned int nEntries12 = curTree->GetEntries();
0567     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
0568       curTree->GetEntry(iEntry);
0569       for (int ii = 0; ii < 12; ++ii) {
0570         surfDeform_[ii] = -1.0;
0571       }
0572       for (int npar = 0; npar < int(inputDpar2.size()); ++npar) {
0573         surfDeform_[npar] = inputDpar2.at(npar);
0574         std::stringstream histname0;
0575         histname0 << "SurfDeform_Par_" << npar;
0576         if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0577           m_h1_[histname0.str()]->Fill(surfDeform_[npar]);
0578         if (inputSubdetid2 == 1) {
0579           std::stringstream histname1;
0580           histname1 << "SurfDeform_PixelBarrel_Par_" << npar;
0581           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0582             m_h1_[histname1.str()]->Fill(surfDeform_[npar]);
0583         }
0584         if (inputSubdetid2 == 2) {
0585           std::stringstream histname2;
0586           histname2 << "SurfDeform_PixelEndcap_Par_" << npar;
0587           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0588             m_h1_[histname2.str()]->Fill(surfDeform_[npar]);
0589         }
0590         if (phase_ < SiPixelPI::phase::two) {
0591           if (phase_ == SiPixelPI::phase::zero) {
0592             (m_vtkmap_.at(npar)).fill_current_val(inputRawid2, surfDeform_[npar]);
0593           } else {
0594             // old-style tracker map does not support Phase1 Pixel
0595             if (inputSubdetid2 != 1 && inputSubdetid2 != 2) {
0596               (m_vtkmap_.at(npar)).fill_current_val(inputRawid2, surfDeform_[npar]);
0597             }
0598           }
0599         }
0600       }
0601     }
0602 
0603   } else if (inputFilename1_ != "IDEAL" && inputFilename2_ == "IDEAL") {
0604     int inputRawid1;
0605     int inputSubdetid1;
0606     int inputDtype1;
0607     std::vector<double> inputDpar1;
0608     std::vector<double>* p_inputDpar1 = &inputDpar1;
0609 
0610     TTree* refTree = (TTree*)inputRootFile1_->Get(inputTreenameDeform_.c_str());
0611     refTree->SetBranchAddress("irawid", &inputRawid1);
0612     refTree->SetBranchAddress("subdetid", &inputSubdetid1);
0613     refTree->SetBranchAddress("dtype", &inputDtype1);
0614     refTree->SetBranchAddress("dpar", &p_inputDpar1);
0615     unsigned int nEntries11 = refTree->GetEntries();
0616 
0617     for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
0618       refTree->GetEntry(iEntry);
0619       for (int ii = 0; ii < 12; ++ii) {
0620         surfDeform_[ii] = -1.0;
0621       }
0622       for (int npar = 0; npar < int(inputDpar1.size()); ++npar) {
0623         surfDeform_[npar] = -inputDpar1.at(npar);
0624         std::stringstream histname0;
0625         histname0 << "SurfDeform_Par_" << npar;
0626         if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0627           m_h1_[histname0.str()]->Fill(surfDeform_[npar]);
0628         if (inputSubdetid1 == 1) {
0629           std::stringstream histname1;
0630           histname1 << "SurfDeform_PixelBarrel_Par_" << npar;
0631           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0632             m_h1_[histname1.str()]->Fill(surfDeform_[npar]);
0633         }
0634         if (inputSubdetid1 == 2) {
0635           std::stringstream histname2;
0636           histname2 << "SurfDeform_PixelEndcap_Par_" << npar;
0637           if (TMath::Abs(surfDeform_[npar]) > (m_rangeHigh_ - m_rangeLow_) / (10. * m_nBins_))
0638             m_h1_[histname2.str()]->Fill(surfDeform_[npar]);
0639         }
0640         if (phase_ < SiPixelPI::phase::two) {
0641           if (phase_ == SiPixelPI::phase::zero) {
0642             (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0643           } else {
0644             // old-style tracker map does not support Phase1 Pixel
0645             if (inputSubdetid1 != 1 && inputSubdetid1 != 2) {
0646               (m_vtkmap_.at(npar)).fill_current_val(inputRawid1, surfDeform_[npar]);
0647             }
0648           }
0649         }
0650       }
0651     }
0652 
0653   } else if (inputFilename1_ == "IDEAL" && inputFilename2_ == "IDEAL") {
0654     edm::LogInfo("TrackerGeometryCompare") << ">>>> Comparing IDEAL with IDEAL: nothing to do! <<<<\n";
0655   }
0656   return;
0657 }
0658 
0659 void TrackerGeometryCompare::compareGeometries(Alignable* refAli,
0660                                                Alignable* curAli,
0661                                                const TrackerTopology* tTopo,
0662                                                const edm::EventSetup& iSetup) {
0663   using namespace align;
0664 
0665   const auto& refComp = refAli->components();
0666   const auto& curComp = curAli->components();
0667 
0668   unsigned int nComp = refComp.size();
0669   //only perform for designate levels
0670   bool useLevel = false;
0671   for (unsigned int i = 0; i < m_theLevels.size(); ++i) {
0672     if (refAli->alignableObjectId() == m_theLevels[i])
0673       useLevel = true;
0674   }
0675 
0676   //another added level for difference between det and detunit
0677   //if ((refAli->alignableObjectId()==2)&&(nComp == 1)) useLevel = false;
0678 
0679   //coordinate matching, etc etc
0680   if (useLevel) {
0681     DetId detid(refAli->id());
0682 
0683     CLHEP::Hep3Vector Rtotal, Wtotal, lRtotal, lWtotal;
0684     Rtotal.set(0., 0., 0.);
0685     Wtotal.set(0., 0., 0.);
0686     lRtotal.set(0., 0., 0.);
0687     lWtotal.set(0., 0., 0.);
0688 
0689     bool converged = false;
0690 
0691     AlgebraicVector diff, check;
0692 
0693     for (int i = 0; i < 100; i++) {
0694       // Get differences between alignments for rotations and translations
0695       // both local and global
0696       diff = align::diffAlignables(refAli, curAli, weightBy_, weightById_, weightByIdVector_);
0697 
0698       // 'diffAlignables' returns 'refAli - curAli' for translations and 'curAli - refAli' for rotations.
0699       // The plan is to unify this at some point, but a simple change of the sign for one of them was postponed
0700       // to do some further checks to understand the rotations better
0701       //Updated July 2018: as requested the sign in the translations has been changed to match the one in rotations. A test was done to change the diffAlignables function and solve the issue there, but proved quite time consuming. To unify the sign convention in the least amount of time the choice was made to change the sign here.
0702       CLHEP::Hep3Vector dR(-diff[0], -diff[1], -diff[2]);
0703       CLHEP::Hep3Vector dW(diff[3], diff[4], diff[5]);
0704       CLHEP::Hep3Vector dRLocal(-diff[6], -diff[7], -diff[8]);
0705       CLHEP::Hep3Vector dWLocal(diff[9], diff[10], diff[11]);
0706 
0707       // Translations
0708       Rtotal += dR;
0709       lRtotal += dRLocal;
0710 
0711       //Rotations
0712       CLHEP::HepRotation rot(Wtotal.unit(), Wtotal.mag());
0713       CLHEP::HepRotation drot(dW.unit(), dW.mag());
0714       rot *= drot;
0715       Wtotal.set(rot.axis().x() * rot.delta(), rot.axis().y() * rot.delta(), rot.axis().z() * rot.delta());
0716 
0717       CLHEP::HepRotation rotLocal(lWtotal.unit(), lWtotal.mag());
0718       CLHEP::HepRotation drotLocal(dWLocal.unit(), dWLocal.mag());
0719       rotLocal *= drotLocal;
0720       lWtotal.set(rotLocal.axis().x() * rotLocal.delta(),
0721                   rotLocal.axis().y() * rotLocal.delta(),
0722                   rotLocal.axis().z() * rotLocal.delta());
0723 
0724       // Move current alignable by shift and check if difference
0725       // is smaller than tolerance value
0726       // if true, break the loop
0727       align::moveAlignable(curAli, diff);
0728       float tolerance = 1e-7;
0729       check = align::diffAlignables(refAli, curAli, weightBy_, weightById_, weightByIdVector_);
0730       align::GlobalVector checkR(check[0], check[1], check[2]);
0731       align::GlobalVector checkW(check[3], check[4], check[5]);
0732       if ((checkR.mag() < tolerance) && (checkW.mag() < tolerance)) {
0733         converged = true;
0734         break;
0735       }
0736     }
0737 
0738     // give an exception if difference has not fallen below tolerance level
0739     // i.e. method has not converged
0740     if (!converged) {
0741       edm::LogInfo("TrackerGeometryCompare")
0742           << "Tolerance Exceeded!(alObjId: " << refAli->alignableObjectId()
0743           << ", rawId: " << refAli->geomDetId().rawId() << ", subdetId: " << detid.subdetId() << "): " << diff << check;
0744       throw cms::Exception("Tolerance in TrackerGeometryCompare exceeded");
0745     }
0746     AlgebraicVector TRtot(12);
0747     // global
0748     TRtot(1) = Rtotal.x();
0749     TRtot(2) = Rtotal.y();
0750     TRtot(3) = Rtotal.z();
0751     TRtot(4) = Wtotal.x();
0752     TRtot(5) = Wtotal.y();
0753     TRtot(6) = Wtotal.z();
0754     // local
0755     TRtot(7) = lRtotal.x();
0756     TRtot(8) = lRtotal.y();
0757     TRtot(9) = lRtotal.z();
0758     TRtot(10) = lWtotal.x();
0759     TRtot(11) = lWtotal.y();
0760     TRtot(12) = lWtotal.z();
0761 
0762     fillTree(refAli, TRtot, tTopo, iSetup);
0763   }
0764 
0765   // another added level for difference between det and detunit
0766   for (unsigned int i = 0; i < nComp; ++i)
0767     compareGeometries(refComp[i], curComp[i], tTopo, iSetup);
0768 }
0769 
0770 void TrackerGeometryCompare::setCommonTrackerSystem() {
0771   edm::LogInfo("TrackerGeometryCompare") << "Setting Common Tracker System....";
0772 
0773   // DM_534??AlignableObjectId dummy;
0774   // DM_534??_commonTrackerLevel = dummy.nameToType(_setCommonTrackerSystem);
0775   commonTrackerLevel_ = currentTracker->objectIdProvider().stringToId(setCommonTrackerSystem_);
0776 
0777   diffCommonTrackerSystem(referenceTracker, currentTracker);
0778 
0779   align::EulerAngles dOmega(3);
0780   dOmega[0] = TrackerCommonR_.x();
0781   dOmega[1] = TrackerCommonR_.y();
0782   dOmega[2] = TrackerCommonR_.z();
0783   align::RotationType rot = align::toMatrix(dOmega);
0784   align::GlobalVector theR = TrackerCommonT_;
0785 
0786   edm::LogInfo("TrackerGeometryCompare") << "what we get from overlaying the pixels..." << theR << ", " << rot;
0787 
0788   //transform to the Tracker System
0789   align::PositionType trackerCM = currentTracker->globalPosition();
0790   align::GlobalVector cmDiff(
0791       trackerCM.x() - TrackerCommonCM_.x(), trackerCM.y() - TrackerCommonCM_.y(), trackerCM.z() - TrackerCommonCM_.z());
0792 
0793   edm::LogInfo("TrackerGeometryCompare") << "Pixel CM: " << TrackerCommonCM_ << ", tracker CM: " << trackerCM;
0794 
0795   //adjust translational difference factoring in different rotational CM
0796   //needed because rotateInGlobalFrame is about CM of alignable, not Tracker
0797   const align::GlobalVector::BasicVectorType& lpvgf = cmDiff.basicVector();
0798   align::GlobalVector moveV(rot.multiplyInverse(lpvgf) - lpvgf);
0799   align::GlobalVector theRprime(theR + moveV);
0800 
0801   AlgebraicVector TrackerCommonTR(6);
0802   TrackerCommonTR(1) = theRprime.x();
0803   TrackerCommonTR(2) = theRprime.y();
0804   TrackerCommonTR(3) = theRprime.z();
0805   TrackerCommonTR(4) = TrackerCommonR_.x();
0806   TrackerCommonTR(5) = TrackerCommonR_.y();
0807   TrackerCommonTR(6) = TrackerCommonR_.z();
0808 
0809   edm::LogInfo("TrackerGeometryCompare") << "and after the transformation: " << TrackerCommonTR;
0810   align::moveAlignable(currentTracker, TrackerCommonTR);
0811 }
0812 
0813 void TrackerGeometryCompare::diffCommonTrackerSystem(Alignable* refAli, Alignable* curAli) {
0814   const auto& refComp = refAli->components();
0815   const auto& curComp = curAli->components();
0816   unsigned int nComp = refComp.size();
0817   //only perform for designate levels
0818   bool useLevel = false;
0819   if (refAli->alignableObjectId() == commonTrackerLevel_)
0820     useLevel = true;
0821 
0822   //useLevel = false;
0823   if (useLevel) {
0824     CLHEP::Hep3Vector Rtotal, Wtotal;
0825     Rtotal.set(0., 0., 0.);
0826     Wtotal.set(0., 0., 0.);
0827 
0828     AlgebraicVector diff = align::diffAlignables(refAli, curAli, weightBy_, weightById_, weightByIdVector_);
0829     CLHEP::Hep3Vector dR(diff[0], diff[1], diff[2]);
0830     Rtotal += dR;
0831     CLHEP::Hep3Vector dW(diff[3], diff[4], diff[5]);
0832     CLHEP::HepRotation rot(Wtotal.unit(), Wtotal.mag());
0833     CLHEP::HepRotation drot(dW.unit(), dW.mag());
0834     rot *= drot;
0835     Wtotal.set(rot.axis().x() * rot.delta(), rot.axis().y() * rot.delta(), rot.axis().z() * rot.delta());
0836 
0837     TrackerCommonT_ = align::GlobalVector(Rtotal.x(), Rtotal.y(), Rtotal.z());
0838     TrackerCommonR_ = align::GlobalVector(Wtotal.x(), Wtotal.y(), Wtotal.z());
0839     TrackerCommonCM_ = curAli->globalPosition();
0840 
0841   } else {
0842     for (unsigned int i = 0; i < nComp; ++i)
0843       diffCommonTrackerSystem(refComp[i], curComp[i]);
0844   }
0845 }
0846 
0847 void TrackerGeometryCompare::fillTree(Alignable* refAli,
0848                                       const AlgebraicVector& diff,
0849                                       const TrackerTopology* tTopo,
0850                                       const edm::EventSetup& iSetup) {
0851   //Get bad modules
0852   const SiPixelQuality* SiPixelModules = &iSetup.getData(pixQualityToken_);
0853   const SiStripQuality* SiStripModules = &iSetup.getData(stripQualityToken_);
0854 
0855   id_ = refAli->id();
0856 
0857   badModuleQuality_ = 0;
0858   //check if module has a bad quality tag
0859   if (SiPixelModules->IsModuleBad(id_)) {
0860     badModuleQuality_ = 1;
0861   }
0862   if (SiStripModules->IsModuleBad(id_)) {
0863     badModuleQuality_ = 1;
0864   }
0865 
0866   //check if module is in a given list of bad/untouched etc. modules
0867   inModuleList_ = 0;
0868   for (unsigned int i = 0; i < moduleList_.size(); i++) {
0869     if (moduleList_[i] == id_) {
0870       inModuleList_ = 1;
0871       break;
0872     }
0873   }
0874 
0875   level_ = refAli->alignableObjectId();
0876   //need if ali has no mother
0877   if (refAli->mother()) {
0878     mid_ = refAli->mother()->geomDetId().rawId();
0879     mlevel_ = refAli->mother()->alignableObjectId();
0880   } else {
0881     mid_ = -1;
0882     mlevel_ = -1;
0883   }
0884   DetId detid(id_);
0885   sublevel_ = detid.subdetId();
0886   fillIdentifiers(sublevel_, id_, tTopo);
0887   xVal_ = refAli->globalPosition().x();
0888   yVal_ = refAli->globalPosition().y();
0889   zVal_ = refAli->globalPosition().z();
0890   align::GlobalVector vec(xVal_, yVal_, zVal_);
0891   rVal_ = vec.perp();
0892   phiVal_ = vec.phi();
0893   etaVal_ = vec.eta();
0894   align::RotationType rot = refAli->globalRotation();
0895   align::EulerAngles eulerAngles = align::toAngles(rot);
0896   alphaVal_ = eulerAngles[0];
0897   betaVal_ = eulerAngles[1];
0898   gammaVal_ = eulerAngles[2];
0899   // global
0900   dxVal_ = diff[0];
0901   dyVal_ = diff[1];
0902   dzVal_ = diff[2];
0903   // local
0904   duVal_ = diff[6];
0905   dvVal_ = diff[7];
0906   dwVal_ = diff[8];
0907   //...TODO...
0908   align::GlobalVector g(dxVal_, dyVal_, dzVal_);
0909   //getting dR and dPhi
0910   align::GlobalVector vRef(xVal_, yVal_, zVal_);
0911   align::GlobalVector vCur(xVal_ + dxVal_, yVal_ + dyVal_, zVal_ + dzVal_);
0912   drVal_ = vCur.perp() - vRef.perp();
0913   dphiVal_ = vCur.phi() - vRef.phi();
0914   // global
0915   dalphaVal_ = diff[3];
0916   dbetaVal_ = diff[4];
0917   dgammaVal_ = diff[5];
0918   // local
0919   daVal_ = diff[9];
0920   dbVal_ = diff[10];
0921   dgVal_ = diff[11];
0922   //detIdFlag
0923   if (refAli->alignableObjectId() == align::AlignableDetUnit) {
0924     if (detIdFlag_) {
0925       if ((passIdCut(refAli->id())) || (passIdCut(refAli->mother()->id()))) {
0926         useDetId_ = 1;
0927       } else {
0928         useDetId_ = 0;
0929       }
0930     }
0931   }
0932   // det module dimension
0933   if (refAli->alignableObjectId() == align::AlignableDetUnit) {
0934     if (refAli->mother()->alignableObjectId() != align::AlignableDet)
0935       detDim_ = 1;
0936     else if (refAli->mother()->alignableObjectId() == align::AlignableDet)
0937       detDim_ = 2;
0938   } else
0939     detDim_ = 0;
0940 
0941   surWidth_ = refAli->surface().width();
0942   surLength_ = refAli->surface().length();
0943   align::RotationType rt = refAli->globalRotation();
0944   surRot_[0] = rt.xx();
0945   surRot_[1] = rt.xy();
0946   surRot_[2] = rt.xz();
0947   surRot_[3] = rt.yx();
0948   surRot_[4] = rt.yy();
0949   surRot_[5] = rt.yz();
0950   surRot_[6] = rt.zx();
0951   surRot_[7] = rt.zy();
0952   surRot_[8] = rt.zz();
0953 
0954   //Fill
0955   alignTree_->Fill();
0956 }
0957 
0958 void TrackerGeometryCompare::surveyToTracker(AlignableTracker* ali,
0959                                              Alignments* alignVals,
0960                                              AlignmentErrorsExtended* alignErrors) {
0961   //getting the right alignables for the alignment record
0962   auto detPB = ali->pixelHalfBarrelGeomDets();
0963   auto detPEC = ali->pixelEndcapGeomDets();
0964   auto detTIB = ali->innerBarrelGeomDets();
0965   auto detTID = ali->TIDGeomDets();
0966   auto detTOB = ali->outerBarrelGeomDets();
0967   auto detTEC = ali->endcapGeomDets();
0968 
0969   align::Alignables allGeomDets;
0970   std::copy(detPB.begin(), detPB.end(), std::back_inserter(allGeomDets));
0971   std::copy(detPEC.begin(), detPEC.end(), std::back_inserter(allGeomDets));
0972   std::copy(detTIB.begin(), detTIB.end(), std::back_inserter(allGeomDets));
0973   std::copy(detTID.begin(), detTID.end(), std::back_inserter(allGeomDets));
0974   std::copy(detTOB.begin(), detTOB.end(), std::back_inserter(allGeomDets));
0975   std::copy(detTEC.begin(), detTEC.end(), std::back_inserter(allGeomDets));
0976 
0977   align::Alignables rcdAlis;
0978   for (const auto& i : allGeomDets) {
0979     if (i->components().size() == 1) {
0980       rcdAlis.push_back(i);
0981     } else if (i->components().size() > 1) {
0982       rcdAlis.push_back(i);
0983       const auto& comp = i->components();
0984       for (const auto& j : comp)
0985         rcdAlis.push_back(j);
0986     }
0987   }
0988 
0989   //turning them into alignments
0990   for (const auto& k : rcdAlis) {
0991     const SurveyDet* surveyInfo = k->survey();
0992     const align::PositionType& pos(surveyInfo->position());
0993     align::RotationType rot(surveyInfo->rotation());
0994     CLHEP::Hep3Vector clhepVector(pos.x(), pos.y(), pos.z());
0995     CLHEP::HepRotation clhepRotation(
0996         CLHEP::HepRep3x3(rot.xx(), rot.xy(), rot.xz(), rot.yx(), rot.yy(), rot.yz(), rot.zx(), rot.zy(), rot.zz()));
0997     AlignTransform transform(clhepVector, clhepRotation, k->id());
0998     AlignTransformErrorExtended transformError(CLHEP::HepSymMatrix(3, 1), k->id());
0999     alignVals->m_align.push_back(transform);
1000     alignErrors->m_alignError.push_back(transformError);
1001   }
1002 
1003   // to get the right order, sort by rawId
1004   std::sort(alignVals->m_align.begin(), alignVals->m_align.end());
1005   std::sort(alignErrors->m_alignError.begin(), alignErrors->m_alignError.end());
1006 }
1007 
1008 void TrackerGeometryCompare::addSurveyInfo(Alignable* ali) {
1009   const auto& comp = ali->components();
1010 
1011   unsigned int nComp = comp.size();
1012 
1013   for (unsigned int i = 0; i < nComp; ++i)
1014     addSurveyInfo(comp[i]);
1015 
1016   const SurveyError& error = theSurveyErrors->m_surveyErrors[theSurveyIndex];
1017 
1018   if (ali->geomDetId().rawId() != error.rawId() || ali->alignableObjectId() != error.structureType()) {
1019     throw cms::Exception("DatabaseError") << "Error reading survey info from DB. Mismatched id!";
1020   }
1021 
1022   const CLHEP::Hep3Vector& pos = theSurveyValues->m_align[theSurveyIndex].translation();
1023   const CLHEP::HepRotation& rot = theSurveyValues->m_align[theSurveyIndex].rotation();
1024 
1025   AlignableSurface surf(
1026       align::PositionType(pos.x(), pos.y(), pos.z()),
1027       align::RotationType(rot.xx(), rot.xy(), rot.xz(), rot.yx(), rot.yy(), rot.yz(), rot.zx(), rot.zy(), rot.zz()));
1028 
1029   surf.setWidth(ali->surface().width());
1030   surf.setLength(ali->surface().length());
1031 
1032   ali->setSurvey(new SurveyDet(surf, error.matrix()));
1033 
1034   ++theSurveyIndex;
1035 }
1036 
1037 bool TrackerGeometryCompare::passIdCut(uint32_t id) {
1038   bool pass = false;
1039   int nEntries = detIdFlagVector_.size();
1040 
1041   for (int i = 0; i < nEntries; i++) {
1042     if (detIdFlagVector_[i] == id)
1043       pass = true;
1044   }
1045 
1046   return pass;
1047 }
1048 
1049 void TrackerGeometryCompare::fillIdentifiers(int subdetlevel, int rawid, const TrackerTopology* tTopo) {
1050   switch (subdetlevel) {
1051     case 1: {
1052       identifiers_[0] = tTopo->pxbModule(rawid);
1053       identifiers_[1] = tTopo->pxbLadder(rawid);
1054       identifiers_[2] = tTopo->pxbLayer(rawid);
1055       identifiers_[3] = 999;
1056       identifiers_[4] = 999;
1057       identifiers_[5] = 999;
1058       break;
1059     }
1060     case 2: {
1061       identifiers_[0] = tTopo->pxfModule(rawid);
1062       identifiers_[1] = tTopo->pxfPanel(rawid);
1063       identifiers_[2] = tTopo->pxfBlade(rawid);
1064       identifiers_[3] = tTopo->pxfDisk(rawid);
1065       identifiers_[4] = tTopo->pxfSide(rawid);
1066       identifiers_[5] = 999;
1067       break;
1068     }
1069     case 3: {
1070       identifiers_[0] = tTopo->tibModule(rawid);
1071       identifiers_[1] = tTopo->tibStringInfo(rawid)[0];
1072       identifiers_[2] = tTopo->tibStringInfo(rawid)[1];
1073       identifiers_[3] = tTopo->tibStringInfo(rawid)[2];
1074       identifiers_[4] = tTopo->tibLayer(rawid);
1075       identifiers_[5] = 999;
1076       break;
1077     }
1078     case 4: {
1079       identifiers_[0] = tTopo->tidModuleInfo(rawid)[0];
1080       identifiers_[1] = tTopo->tidModuleInfo(rawid)[1];
1081       identifiers_[2] = tTopo->tidRing(rawid);
1082       identifiers_[3] = tTopo->tidWheel(rawid);
1083       identifiers_[4] = tTopo->tidSide(rawid);
1084       identifiers_[5] = 999;
1085       break;
1086     }
1087     case 5: {
1088       identifiers_[0] = tTopo->tobModule(rawid);
1089       identifiers_[1] = tTopo->tobRodInfo(rawid)[0];
1090       identifiers_[2] = tTopo->tobRodInfo(rawid)[1];
1091       identifiers_[3] = tTopo->tobLayer(rawid);
1092       identifiers_[4] = 999;
1093       identifiers_[5] = 999;
1094       break;
1095     }
1096     case 6: {
1097       identifiers_[0] = tTopo->tecModule(rawid);
1098       identifiers_[1] = tTopo->tecRing(rawid);
1099       identifiers_[2] = tTopo->tecPetalInfo(rawid)[0];
1100       identifiers_[3] = tTopo->tecPetalInfo(rawid)[1];
1101       identifiers_[4] = tTopo->tecWheel(rawid);
1102       identifiers_[5] = tTopo->tecSide(rawid);
1103       break;
1104     }
1105     default: {
1106       edm::LogInfo("TrackerGeometryCompare") << "Error: bad subdetid!!";
1107       break;
1108     }
1109   }
1110 }
1111 
1112 DEFINE_FWK_MODULE(TrackerGeometryCompare);