OfflineValidation/plugins/TrackerGeometryCompare.cc

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